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INTRODUCTION 

BACKGROUND 

Until  recently,  the  aerodynamic  analysis  of  helicopter  configurations  has 
been  largely  neglected,  with  the  exception  being  the  analysis  of  rotor 
systems.  Even  in  this  area,  a particular  rotor  airfoil  cross-section  is 
the  result  of  many  hours  of  exhaustive  wind  tunnel  testing,  usually  with 
only  a minimal  amount  of  theoretical  support.  This  pattern  in  changing, 
however,  due  to  a demand  for  greater  vehicle  performance,  either  in  speed, 
lifting  capacity  or  endurance.  Another  consideration,  the  more  efficient 
use  of  energy,  has  further  altered  the  traditional  view  of  helicopter  aero- 
dynamicists.  Drag  has  become  an  important  issue,  and  work  is  now  under  way 
on  several  fronts  in  an  attempt  to  understand  and  subsequently  reduce  its 
effect  through  improved  helicopter  design. 

The  analysis  of  current  helicopter  fuselage  designs  is  basically  a problem 
in  the  calculation  of  three-dimensional  bluff  body  aerodynamics.  Flow 
separation  results  in  a flow  field  problem  which,  while  not  necessarily 
intractable,  remains  one  of  the  most  challenging  of  fluid  dynamic  phenom- 
ena. Recourse  to  the  wind  tunnel  is  only  moderately  successful  because  of 
tunnel  blockage  and  Reynolds  number  effects.  Theoretical  methods  fail 
through  inability  to  model  separation  effects,  boundary  layer  displacement 
effects,  or  the  interaction  of  the  rotor  wake  with  the  fuselage  pressure 
field.  While  it  appears  that  the  problems  inherent  in  wind  tunnel  testing 
will  remain  with  us,  improvements  in  computer  simulation  may  ultimately 
lead  to  analytical  methods  which  will  supplant  the  wind  tunnel.  In  fact, 
a recent  paper  by  Chapman,  Mark  and  Pirtle  (Reference  1)  of  the  NASA-Ames 
Research  Center,  projects  that  this  will  occur  within  a decade  if  the  cur- 
rent rate  of  advancement  in  computer  technology  continues.  Although  the 
complete  solution  of  the  Navier  Stokes  equations  for  flow  around  an  arbi- 
trary body  remains  beyond  our  present  capability,  much  useful  work  has 
been  done  in  the  development  of  three-dimensional  aerodynamic  analysis 
tools.  The  ability  to  model  separation  would  greatly  enhance  the  useful- 
ness of  these  methods. 

Several  approaches  have  been  taken  to  model  bluff  body  separation,  i.e., 
by  filling  in  the  separated  region  with  sources  or  a solid  body  fairing 
(Reference  2) . These  methods  cannot  predict  the  separation  pressure 
levels  directly,  and  consequently  rely  on  some  criterion  to  specify  the 
pressure  level.  They  do,  however,  generally  improve  the  predicted  pres- 
sure distributions  upstream  of  separation.  Improved  performance  predic- 
tion requires  that  the  separated  region  pressure  level  be  predicted  as 
part  of  the  overall  analysis  without  recourse  to  empiricism.  A separa- 
tion model  based  on  a vorticity  shear  layer  approach  which  meets  the 
above  requirement  has  been  developed  and  incorporated  into  a general 


three-dimensional  potential  flow  program.  The  separation  model  and  the  re- 
sulting analysis  procedure  are  described  in  the  following  sections. 

GENERAL  APPROACH 

An  analysis  procedure  capable  of  predicting  the  three-dimensional  aerodynam- 
ic characteristics  of  helicopter  fuselage  configurations  having  separated 
flow  regions  was  the  objective  of  the  proposed  study.  In  accomplishing 
this  objective,  two  major  tasks  were  completed: 

(1)  The  modification  of  Computer  program  WBAERO  in  order  to 
include  viscous  effects  in  the  calculation  of  the  potential 
flow.  In  this  case,  viscous  effects  include  the  displacement 
of  external  streamlines  both  by  the  upstream  attached  boundary 
layer  and  by  the  separated  flow  region.  In  many  cases,  a very 
thick  but  attached  boundary  layer  will  displace  the  external 
flow  nearly  as  much  as  would  a separation  region. 

(2)  The  development  and  verification  of  the  vorticity  separation 
model  described  in  this  report. 

The  general  analysis  procedure  consists  of  a potential  flow  calculation,  a 
boundary  layer  analysis,  and  a model  of  the  separated  flow.  The  analysis 
begins  with  the  creation  of  planar  surface  panels  from  fuselage  cross-sec- 
tion data  input  by  the  user.  The  inviscid  flow  field  about  this  configura- 
tion is  determined  by  the  potential  flow  program.  A series  of  streamlines 
are  then  calculated  over  the  configuration  surface,  providing  information 
required  for  the  boundary  layer  methods.  Laminar,  transition,  and  turbulent 
boundary  layer  characteristics  (including  the  separation  point,  if  separa- 
tion is  present)  are  determined  along  each  streamline,  as  are  the  source 
distributions  representing  the  attached  flow  viscous  effects . The  separa- 
tion points  are  used  to  map  the  separated  flow  region  on  the  configuration 
geometry  from  which  the  vorticity  shear  layer  leaves  the  surface.  A new 
potential  flow  calculation  is  made  with  the  inclusion  of  both  the  separa- 
tion model  and  the  boundary  layer  sources.  The  result  is  a predicted  pres- 
sure field  for  the  entire  body  in  both  attached  and  separated  flow. 
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POTENTIAL  FLOW  METHOD 


The  potential  flow  on  the  body  is  computed  using  the  WBAERO  program  des- 
cribed in  Reference  3.  The  body  surface  is  divided  into  a large  number  of 
planar  panels,  each  of  which  contains  a constant  source  distribution.  In 
addition,  an  internal  vortex  lattice  system  may  be  located  along  the  mean 
chord  plane  of  lifting  surfaces  to  provide  circulation  to  the  flow.  A 
typical  helicopter  fuselage  panel  subdivision  is  shown  in  Figure  1. 

30METRY  DEFINITION 


The  body  is  described  by  a series  of  cross-sections  given  at  selected 
intervals  along  its  length.  The  surface  panels  are  located  between  adja- 
cent sections,  with  the  corner  points  being  defined  by  the  input  cross- 
section  coordinates.  Although  the  standard  WBAERO  geometry  input  has  been 
retained  in  the  current  program,  it  can  easily  be  adapted  to  incorporate 
more  sophisticated  geometry  inputs  such  as  the  Sikorsky  Fuselage  Geometry 
Definition  Program,  Reference  4. 

INVISCID  METHOD 


Analytical  expressions  for  the  perturbation  velocity  field  induced  by  a 
constant  source  distribution  on  an  arbitrary  quadrilateral  panel  are  given 
by  Hess  and  Smith  (Reference  5) . Similarly,  the  velocity  field  induced  by 
the  elements  of  a vortex  lattice  are  given  by  Rubbert  and  Saaris  (Refe- 
rence 6) . The  perturbation  velocities  are  used  to  calculate  the  coeffici- 
ents of  a system  of  linear  equations  relating  the  magnitude  of  the  normal 
velocities  at  the  panel  control  points  to  the  unknown  source  and  vortex 
strengths.  The  source  and  vortex  strengths  which  satisfy  the  boundary 
condition  of  tangential  flow  at  the  control  points  for  a given  Mach  number 
and  angle  of  attack  are  determined  by  solving  this  system  of  equations  by 
an  iterative  procedure.  The  pressure  coefficients  at  panel  control  points 
are  then  calculated  in  terms  of  the  perturbation  velocity  components,  and 
the  forces  and  moments  acting  on  the  configuration  are  obtained  by  numeri- 
cal integration. 

The  perturbation  velocity  components  induced  by  the  sources  and  vortices 
are  described  in  Reference  3,  together  with  the  formation  and  solution 
of  the  boundary  condition  equations,  and  the  procedure  used  to  calculate 
the  pressure  coefficients,  forces,  and  moments  on  the  configuration. 


1 

STREAMLINE  CALCULATION 


The  WBAERO  program  computes  an  approximate  solution  to  the  problem  of 
three-dimensional  potential  flow  over  an  arbitrary  body.  The  body  surface 
is  represented  by  a set  of  plane  quadrilaterals,  and  the  three  components 
of  the  velocity  vector  are  computed  at  the  panel  centroids.  The  stream- 
line program  uses  the  velocity  data  at  four  neighboring  panels  to  provide 
a linear  approximation  to  the  velocity  anywhere  on  a given  panel.  A 
straight-line  approximation  to  the  actual  streamline  passing  through  a 
point  on  the  panel  is  constructed,  and  continued  over  all  panels  upstream 
and  downstream.  This  results  in  relatively  smooth  streamlines  over  the 
body,  and  allows  the  geodesic  curvatures  of  the  streamlines  and  equipoten- 
tial  lines  to  be  computed.  From  this  information,  the  metric  or  effective 
radius  of  the  body  along  the  streamline  is  obtained.  The  method  used  fol- 
lows closely  that  described  in  Reference  7. 

COMPUTATION  OF  THE  PARTIAL  DERIVATIVES 

The  quadrilateral  under  consideration,  and  its  four  neigbors,  is  shown  on 
the  following  sketch.  The  quadrilateral  coordinate  system  is  identified 
by  the  letter  Q,  the  upper  neighbour  by  U,  the  lower  by  L,  the  upstream 
by  T,  and  the  downstream  by  R. 
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The  neighboring  panels  are  just  assumed  to  be  rotated  about  the  common 
edges  so  that  all  panels  lie  in  the  same  plane  as  the  central  quadrila- 
teral. The  program  then  computes  the  velocity  and  position  vectors  of  the 
neighbors  in  terms  of  the  central  panel  coordinate  system.  From  this  con- 
struction, the  following  partial  derivatives  are  obtained. 


(3) 


(4) 


where  VU  is  the  transverse  component  of  velocity  in  the  local  coor- 
dinate system  of  the  upper  neighbor, 

UU  is  the  axial  component  of  velocity  in  the  local  coordinate 
system  of  the  upper  neighbor,  etc., 

DXR  is  the  axial  distance  along  the  surface  between  the  cen- 
troids of  panels  Q and  R, 

DXT  is  the  axial  distance  along  the  surface  bwtween  the  cen- 
troids of  panels  Q and  T, 

DYLJ  is  the  transverse  distance  along  the  surface  between  the 
centroids  of  panels  Q and  U,  and 

DYL  is  the  transverse  distance  along  the  surface  between  the 
centroids  of  panels  Q and  L. 

The  axial  and  transverse  components  of  velocity  on  panel  Q are  then  ap- 
proximately: 

U = Uq  + UjXq  + U2yq  V = Vq  + VlXq  + V2yq  (5) 
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COMPUTATION  OF  POINTS  ON  A STREAMLINE 

The  streamline  computation  begins  at  a given  point  in  a specified  quadri- 
lateral with  the  velocity  vector,  V = (U,V) . A stream  function  is  desired 
which  will  have  a constant  value  of  zero  along  the  streamline.  Since  the 
divergence  of  ^ is  not  zero,  an  additional  function,  p,  must  be  found  such 
that  the  divergence  of  pV  is  zero.  Then  for  a particle  following  a stream- 
line 

= dy/dt  = V = (6) 

dx  dx/dt  U nU 


Thus  a new  vector  field  is  constructed  whose  streamlines  are  identical  to 
those  of  the  velocity  field,  but  whose  divergence  is  zero.  The  function, 
p,  is  given  by  the  series 


U (U  + V ) 
p(x  , y ) = 1 3 L L 

q q U 2 + V 2 

q q 


v (u  + v ) 
q l 2 


U 2 + V 2 

q q 


y + . . . 
q 


The  resulting  stream  function  is 


SF(x  ,y  ) = SF  - V x + U v 

q q 0 q q qf  q 


U V (U  + V ) x 2 
V - -3  3.  .1 L _R_ 

1 u 2 + v 2 2 

q q 


u v (u,  + v ) 1 y 
q q 1 2 -C 


2 2 2 
U + V 

q q 


u (u  + v ) 

q 1 2 


u 2 + V 2 

q q 


Vq  + 


The  constant  value,  SF^,  is  chosen  so  that  the  stream  function  is  zero  at 
the  specified  point.  Since  the  velocity  is  known  only  through  the  linear 
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terms,  it  is  consistent  to  truncate  the  stream  function  after  the  quadratic 
terms . 

Values  of  SF  at  the  four  corner  points  of  the  quadrilateral  are  now  compu- 
ted. By  comparing  the  signs  of  SF  at  adjacent  corners,  the  sides  through 
which  the  streamline  passes  are  determined.  If  the  streamline  passes 
through  a side,  the  value  of  SF  is  computed  at  the  midpoint  of  the  side. 
From  these  three  points,  the  intersection  point  is  computed  from  a three- 
point  interpolation  formula. 


The  parameter,  t(0  < t < 1),  is  defined  so  that 


Here, 


Then 


X = 

X 

1 

+ 

t(x  - X ) 
3 1 

y = 

y 

i 

+ 

t (y  - y ) 

3 1 

x,y 

is 

the 

intersection  point. 

x ,y 
i i 

is 

one 

corner  point. 

x ,y 

3 3 

is 

the 

other  corner  point, 

and 

x ,y 
2 2 

is 

the 

middle  point. 

0 = 

SF (x,y) 

= 2SF  (t  - J*)  (t  - 

1) 

(9) 


+ 2SF  t(t  - V 

3 


(10) 


The  root  of  this  equation  is  chosen  so  that  0 < t < 1.  Note  that  if  there 
had  been  two  roots  between  0 and  1,  then  SF  would  have  the  same  sign  at 
both  corner  points,  and  those  intersection  points  would  be  ignored  com- 
pletely . 

In  general,  two  intersections  will  be  found  for  the  entire  quadrilateral: 
one  where  the  streamline  enters,  and  one  where  it  leaves.  It  is  possible, 
however,  for  there  to  be  four  intersection  points.  The  next  point  on  the 
streamline  is  chosen  from  the  intersection  points  as  follows:  a quantity, 

Q,  is  computed  for  each  intersection  point  by  taking  the  dot  product  of 
the  vector  from  the  starting  point  to  the  intersection  point  with  the 
velocity  vector  at  the  centroid  of  the  quadrilateral,  and  then  multiply- 
ing it  by  the  direction  sign.  (The  direction  sign  is  +1  if  the  stream- 
line is  being  traced  upstream.)  The  intersection  point  with  the  maximum 
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positive  value  of  Q is  chosen.  However,  if  the  largest  Q is  less  than 
or  equal  to  zero,  none  of  the  points  is  acceptable,  and  the  program  search- 
es for  another  quadrilateral  through  which  the  streamline  might  pass. 

Since  each  of  the  points  representing  a streamline  is  located  at  the  bound- 
ary between  two  quadrilaterals,  the  values  of  the  velocity  and  geodesic 
curvatures  at  each  point  can  be  computed  separately  for  the  two  quadrila- 
terals. The  average  of  the  two  values  at  a point  is  taken  as  the  value 
at  that  point. 

The  streamline  is  first  traced  in  the  downstream  direction,  and  then  in 
the  upstream  direction.  After  the  entire  streamline  has  been  obtained, 
the  arc  length  along  the  streamline  is  computed  so  that  it  starts  at  zero 
at  the  upstream  end.  Points  which  are  very  close  together  are  combined. 
Such  points  occur  when  the  streamline  cuts  across  a corner  of  a quadrila- 
*-eral*  distance  between  two  points  is  less  them  one-eighth  the  dis- 

tance between  their  neighbor  in  front  and  their  neighbor  in  back,  they  are 
combined,  and  average  values  of  their  velocity,  position,  and  curvatures 
are  used  for  the  new  point. 

COMPUTATION  OF  THE  GEODESIC  CURVATURES , (K  AND  K ) , AND  THE  METRIC 

- — j 2 

COEFFICIENT,  H 
2 

The  unit  vector  tangent  to  the  streamline  is 

* = IP  .♦  jv  + fa  uu 

(u2  + v2  + w2)** 

are  the  components  of  the  velocity  vector,  and 
are  unit  vectors  in  the  three  coordinate  directions . 


where  U,  V,  W 

-f  -*■ 

i,  j,  k 


The  curvature  of  the  streamline  is  defined  by 

K = dT/di  (12) 

where  Z is  the  arj  length  along  the  streamline.  The  geodesic  curvature  is 
the  component  of  K on  the  body  surface 
■y 

Kg  = K • (T  x N)  (13) 

where  N is  the  unit  normal  vector  to  the  surface. 

For  a plane  quadrilateral  and  with  vectors  written  in  the  quadrilateral 
coordinate  system,  we  have 
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2 2 ^ 
(iU  + jV)/(U  + V P 


(14) 


T = 


d 

d£ 

ius  * 

V )/(U2  -r  V2)*5 

(15) 

■> 

K = 

(tv  - -Jo) 

U(VU  - UV  ) + V(VU  - UV  ) ] /(U2  + V2) 2 

xx  y y J 

(16) 

where 

subscripts 

x and  y 

denote  partial  differentiation 

-►  ■+ 

T x N = 

(-iv  + ju)/(u2  + V2)^ 

Thus , 

K , the  geodesic  curvature  for  the  streamline  is  obtained. 
2 3 

K 

2 

u(uv 

X 

1 ? , 2 

- VU  ) + V(UV  - VU  ) /(U2  + V2) 

x y y 

(17) 

The  geodesic  curvature  for  the  equipotential  lines,  K , is  found  by  simply 

l 

replacing  U by  V and  V by  -U  in  the  equation  for  K 


K 

1 


-V(VU 


UV  ) 
x 


+ U(VU 


UV  ) 
Y 


I /(IT 


_3 

v2) 2 


(18) 


The  metric  coefficient,  H , is  computed  from  the  values  of  K . By  defini- 
tion. 


K 


l 


H H 3 1 


i a 


H is  defined  to  be  unity  along  the  streamline,  so 


H K 
2 l 


3 H 
2 

n 


(19) 


(20) 
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This  equation  is  approximated  by 


or 


H (£  + A£)  - H (£) 
2 2 


H ( £)  K (£) 

2 2 J 


+ 


H (£  + A£)  K (£  + A£) 
2 1 


H (£  + A£) 
2 


H (£) 
2 


2 - A£  K (£) 

I 

2 + A£  K (£  + A£) 

l 


(21) 


where  A£  is  the  chord  length  of  the  segment  of  the  streamline  passing 
through  the  quadrilateral. 

Note  that  H is  determined  only  to  within  a multiplication  constant.  This 
2 

constant  was  arbitrarily  chosen  so  that  is  unity  at  the  specified  start- 
ing point  for  the  streamline.  Thus,  for  an  axisymraetric  body,  will  be 
proportional  to  the  radius,  but  it  will  not,  in  general,  equal  the  radius. 
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BOUNDARY  LAYER  METHODS 


The  boundary  layer  development  on  an  arbitrarily  shaped  three-dimensional 
lifting  configuration  can  be  very  complex.  A thorough  and  exact  calcula- 
tion of  this  development  is  properly  the  domain  of  the  time-dependent  solu- 
tion to  the  general  Navier  Stokes  equations.  Unfortunately,  the  computer 
does  not  yet  exist  which  is  capable  of  handling  such  a problem,  and  even 
if  one  did,  the  cost  in  computer  time  would  be  astronomical.  Less  diffi- 
cult or  costly  are  the  three-dimensional  boundary  layer  programs  now  in 
existence  (References  8,  9 and  10) . These  are,  however,  all  finite  dif- 
ference methods  which  are  complex,  difficult  to  use  (see  Nash,  Reference  8), 
and  generally  very  configuration  dependent.  The  amount  of  computer  time 
required  for  each  calculation  still  prohibits  their  use  in  an  analysis  pro- 
cedure of  the  type  reported  herein.  Having  made  the  above  evaluation,  one 
must  conclude  that  if  the  objective  is  a viscosity-dependent  calculation 
procedure  of  practical  use  to  the  aerodynamicist  for  drag  analysis  and, 
possibly,  for  preliminary  design,  the  method  must  be  relatively  simple  to 
use  and  economic  of  computer  time . This  can  only  be  achieved  if  integral 
boundary  layer  methods  are  used.  In  two  dimensions,  integral  methods  are 
typically  about  100  times  faster  than  finite  difference  methods.  When 
compared  with  three-dimensional  finite  difference  methods,  the  time  dis- 
parity becomes  considerably  greater. 

The  use  of  integral  methods  along  external  streamlines  taking  into  account 
streamline  convergence  or  divergence  has  proven  to  be  a valid  and  reliable 
approach.  They  can,  however,  be  expected  to  break  down  in  regions  of 
large  streamline  convergence  or  divergence  (regions  of  large  crossflow) . 

As  shown  in  the  photographs  of  Reference  2,  this  occurs  usually  in  the  re- 
gion of  separation  where  none  of  the  boundary  layer  methods  (including 
three-dimensional)  can  be  expected  to  be  valid.  It  is  anticipated,  there- 
fore, that  integral  methods  will  suffice  for  most  applications  of  interest 
to  the  aerodynamicist  for  drag  prediction  or  preliminary  design. 

In  those  cases  of  special  interest  to  the  aerodynamicist,  such  as  the  ef- 
fect of  area  suction  for  boundary  layer  control  or  of  roughness  (rivets, 
etc.)  on  drag,  alternate  boundary  layer  calculation  modules  are  available. 
These  methods  are  called  as  needed  into  the  overall  calculation  procedure 
A brief  description  of  all  of  the  boundary  layer  methods  is  given  in  the 
following  paragraphs. 
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LAMINAR  METHOD 


The  integral  method  is  an  adaption  by  Curie  (Reference  11)  of  a method 
developed  by  Thwaites  (Reference  12).  In  Thwaites'  method,  the  momentum 
integral  equation 

de/dx  = Cf/2  - (H  + 2)(0/Xj(dU/dx)  - (6/rXdr/dx)  (22) 


is  written  in  the  form 

Cd/dxXK/U)  = L/U  -(K/r) (dr/dx) / (dU/dx)  (23) 


where 

K = 02/V (dU/dx) 

L = U - K(H  + 2)  } 
l = (0/U)  OU/3y)y=0 


(24) 


1 dr 
r dx 


a measure  of  the  streamline  divergence  or  convergence 


In  Thwaites'  method,  the  divergence  term  was  not  considered,  although  the 
method  was  later  extended  to  include  this  term  by  Rott  and  Crabtree  (Refe- 
rence 13)  . 

Thwaites  used  exact  solutions  to  a variety  of  laminar  flows  to  determine 
the  relationship  between  L and  K, 

L = 0.45  - 6K  (25) 


Curie  has  pointed  out  that  Equation  (25)  is  not  adequate  in  flows  approach- 
ing separation,  and  he  has  suggested  an  extension  or  correction  giving 

L = 0.45  - 6K  + g(K,p)  (26) 

The  parameter,  p,  is  a function  of  both  the  pressure  gradient  and  the  cur- 
vature or  second  derivative  of  velocity. 
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M = K2U(d2U/dx2)/(dU/dx)2  (27) 

Curie  rewrote  Equation  (26)  in  the  form 

L = Fq(K)  - MGq(K)  (28) 


where  Fq  and  Gq  are  universal  functions  determined  from  a series  of  exact 

solutions  to  laminar  flows  in  the  same  way  as  they  are  in  Equation  (25)  . 
After  substitution  of  Equation  (26)  into  Equation  (23),  and  with  subsequent 
integration,  the  result  can  be  rearranged  in  the  form 

x 

02  = 0 . 45v/ (U6r2 ) f r2 (1  + 2.22g)USdx  (29) 

o 

This  equation  is  conveniently  solved  by  iteration,  g initially  equal  to 
zero.  With  values  of  K and  p determined  in  the  first  iteration,  a second 
iteration  is  carried  out  using  Equation  (28).  At  each  step  in  the  calcula- 
tion, the  local  skin  friction  coefficient,  Cf  , and  the  shape  factor,  H, 

can  be  calculated  using  Equation  (24)  . The  local  skin  friction  coefficient 
has  been  defined  as 

Cf  = (U/p0U)£  (30) 

where  £ in  Equation  (24) is  determined  in  a similar  manner  to  L from  a 
series  of  known  solutions  to  give 

l2  = F1(K)  - G^ (K)  (31) 

The  functions  FQ,  F^,  and  G^  are  tabulated  in  the  computer  program. 

Calculations  begin  at  the  stagnation  point,  with  the  initial  momentum 
thickness,  9,  given  as  a function  of  K.  For  bluff  bodies,  K takes  an  ini- 
tial value  (Kq  = 0.0604)  at  the  stagnation  point,  from  which  the  initial 

momentum  thickness,  0 , is 

o 

0 = (.0604  /dU/dx) 1/2  (32) 

o 

The  calculation  proceeds  either  to  laminar  separation  or  to  the  end  of  the 
airfoil  , whichever  occurs  first.  The  calculated  boundary  layer  develop- 
ment is  then  interrogated  to  determine  if  transition,  laminar  separation 
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or  forced  transition  (boundary  layer  tripping)  has  taken  place . If  any  of 
these  phenomena  have  occurred,  the  downstream  flow  is  assumed  to  be  turbu- 
lent. 

Boundary  Layer  Transition  and  Laminar  Separation 

Boundary  layer  transition  is  a very  complex  phenomenon,  and  to  date  no  com- 
pletely reliable  general  theoretical  method  has  been  developed  for  its  pre- 
diction. Reynolds  number  is  a controlling  parameter,  but  it  has  been  shown 
that  the  Reynolds  number  at  transition  can  be  increased  a considerable 
amount  by  careful  elimination  of  disturbances.  At  very  low  Reynolds  num- 
bers, laminar  boundary  layers  are  stable  to  small  disturbances.  However, 
at  higher  Reynolds  numbers,  the  boundary  layer  is  unstable,  and  small  dis- 
turbances can  be  amplified.  Amplification  of  these  disturbances  causes 
the  flow  to  become  turbulent.  The  point  at  which  flow  breakdown  oucurs 
depends  on  the  strength  and  dominant  frequency  of  the  initial  disturbance. 
Disturbances  may  be  due  to  freestream  turbulence,  surface  roughness,  noise 
or  vibration  of  the  surface.  As  there  is  no  detailed  analysis  of  the  tran- 
sition process,  transition  prediction  is  accomplished  by  means  of  empiri- 
cal correlations.  Granville  (Reference  14  ) has  developed  a procedure 
based  on  the  determination  of  the  neutral  stability  point  and  the  tran- 
sition point.  The  neutral  stability  point  is  defined  as  that  point  down- 
stream from  which  small  disturbances  are  amplified  within  the  boundary  layer. 
It  is  this  amplification  of  small  disturbances  that  ultimately  leads  to 
transition.  The  neutral  stability  point  is  reached  when  the  Reynolds 
number,  based  on  the  local  momentum  thickness  and  local  flow  properties, 

attains  some  critical  value , Rn . . S'hlichtinq  and  Ulrich  (Reference  15  ) 

Bins 

have  shown  that  Rglng  can  he  correlated  with  the  local  pressure  gradient 

parameter,  K = (02/v) (dU/ds) . Correlations  by  Smith  (Reference  16)  and 
others  have  been  reduced  to  analytical  form  as  follows. 

Instability  Curves 

K = - 0.4709  + 0.11066  lnRy  - 0.0058591  ln2RQ  (33) 

for  0 < Rq.  < 650 
0ins  — 


and 

K = 0.69412  - 0.23992  lnR0  + 0.0205  in2RQ  (34) 

for  650  < RQ.  < 10,000. 
urns  — 


i 

t 
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If  for  a given  Rg  the  pressure  gradient  parameter,  K,  as  calculated  by 

Equation  (33)  or  (34),  is  greater  than  that  determined  by  the  boundary  layer 
developemnt,  the  flow  has  passed  from  a stable  to  an  unstable  region.  Once 
the  flow  passes  into  the  unstable  region,  the  transition  process  begins; 
and  Granville  has  been  able  to  show  that  a correlation  similar  to  the  in- 


stability process  can  be  used  to  determine  the 
med  an  average  pressure  gradient  parameter,  K, 

s. 


transition  point, 
defined  as 


He  for- 


trans 
ds 


ins 


trans 


(35) 


which  correlated  reasonably  well  with  the  momentum  thickness  Reynolds  num- 
ber at  transition  RQtrans*  This  correlation  is  presented  in  analytical 

form  as  follows. 


Transition  Curves 


K = - 0.0925  + 7.0  x 10~5R„ 

D 

for  0 < Rq,  <750, 

W trans  — 

K = - 0.12571  + 1.14286  x 10_4Ro 

o 

for  750  < Ra  <1,100 

otrans  — 


(36) 


(37) 


and 


K = 1.59381  - 0.45543 


lnRe 


+ 0.032534  ln‘ 


*0 


for  1,100  < Rgtrans  £ 3,000. 


(38) 


When  the  K calculated  by  one  of  the  above  expressions  for  a given  Rg  is 

greater  than  the  value  determined  from  the  boundary  layer  development, 
transition  is  predicted. 

With  transition  predicted,  initial  values  of  the  momentum  thickness,  0, 
and  the  shape  factor,  H,  are  required  to  start  the  turbulent  boundary 
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layer  calculation.  Because  the  boundary  layer  growth  is  continuous,  the 
momentum  thickness  at  transition  is  used  as  the  initial  turbulent  momentum 
thickness.  Since  the  shape  factor  varies  from  values  greater  than  2.0  to 
less  than  1.5  through  the  transition  region,  an  empirical  expression  is 
used  to  determine  the  initial  turbulent  shape  factor.  The  empirical  rela- 


tion between  H and  R, 
ference  17  ) : 


Strans 


was  determined  from  data  obtained  by  Coles  (Re- 


1.4754 


Log10  R0trans 


+ 0.9698 


(39) 


In  many  cases,  the  pressure  gradient  is  of  sufficient  strength  to  separate 
the  laminar  boundary  layer  prior  to  transition.  Except  in  extreme  cases, 
the  boundary  layer  will  then  reattach,  usually  as  a turbulent  boundary 
layer.  Only  recently  have  researchers  been  able  to  analyze  this  phenomenon 
(Reference  18  ) , and  as  yet  the  procedure  is  extremely  complicated  and  cum- 
bersome; consequently,  empirical  relationships  are  still  required.  From 
the  measurements  of  Gaster  (Reference  19  ) and  others,  a correlation  is  formed 
which  is  capable  of  predicting  both  the  occurrence  of  a separation  and 
later  the  reattachment  as  a turbulent  boundary  layer  or  the  catastrophic 
separation.  The  correlation  is  of  the  form 


K = 0.27  - 0.0007575  R0  - 0.000001157  R0 


(40) 


for  R0  > 125 


and 


.09 


(41) 


for  R0  < 125. 


If  K becomes  less  than  -0.09,  separation  occurs,  and  if  R0 
125,  the  boundary  layer  is  not  able  to  reattach.  However  if  R 


is  less  than 

is  greater 


than  125,  the  value  of  K determined  by  the  boundary  layer  development  must 
be  less  than  that  calculated  by  Equation  (40)  before  separation  without 
reattachment  is  predicted.  If  reattachment  is  predicted,  the  turbulent 
boundary  layer  calculation  is  initiated  using  the  momentum  thickness  cal- 
culated at  the  separation  point. 
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TURBULENT  METHOD 


Methods  for  the  calculation  of  turbulent  boundary  layers  in  two  dimensions 
have  been  developed  by  many  investigators.  These  methods  were  reviewed 
at  a conference  held  in  1968  at  Stanford  University  (Reference  20) . 

One  of  the  methods,  an  integral  method  by  Nash  and  Hicks  (Reference  21), 
compared  very  favorably  with  the  more  complex  finite  difference  methods. 
Now,  several  years  later,  the  method  remains  an  excellent  approach  for 
application  to  the  current  problem,  both  in  terms  of  accuracy  and  speed. 

The  Nash-Hicks  method  is  based  on  momentum  and  moment  of  momentum  equa- 
tions coupled  with  a skin  friction  law  derived  from  Coles’  velocity  profile 
family  (Reference  22) . An  additional  equation  is  obtained  by  relating  the 
shear  stress  integral  to  its  equilibrium  value  using  a simple  first-order 
differential  equation.  The  equations  have  been  derived  in  References  23 
and  24,  and  are  repeated  herein  for  completeness. 

A family  of  integral  equations  can  be  derived  taking  higher  moments  of  the 
equation  of  motion.  The  resulting  equations  can  be  expressed  as 


i iil 
r 3y 


6 


3 (ur) 
3* 


dy}  yady 


a + 1 


du 

dx 


o 

0 gives  the  momentum  integral  equation, 

1 gives  the  moment  of  momentum  integral 


and 

equation 


(42) 


The  velocity  distribution  across  the  boundary  layer  can  be  represented  by 
Coles'  velocity  profile  family,  given  by 


u 


u ;ln(y  u ) + Cl  + uQ  , , i 

T 1 •*  T 1 _$_  { 1 - cos  UTy) } 

K V 2 5 


where 


friction  velocity  (T^/P) 2 


u0  = free  parameter  having  units  of  velocity, 

p 

5 = boundary  layer  thickness, 


K = .41,  and 


(43) 
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C = 2.05. 


Substitution  of  this  equation  into  Equation  (42)  results  in  two  equations 
of  the  form: 


du 


dx 


du 

B _& 
dx 


♦ £*  + 
6 d* 


dll 


E dr 


dx 


r dx 


(44) 


A third  equation  of  the  same  form  is  obtained  by  evaluating  Equation  (43) 
at  y = 6 followed  by  the  differential  with  respect  to  x.  The  parameter,  $, 
is  represented  by  the  following  relations: 


(i)  a = 0 

(ii)  a = 1 

(iii)  a = oo 


$ 

$ 


-1 


6 ux 


^1 

iS2 

$ = 0 


( 

J 


dy 


(45) 


The  shear  stress  integral 


o 

J 


P 


dy,  appearing  in  Equation  (45)  was  evalu- 


ated by  Nash  and  Hicks  using  an  equation  of  the  form 
d C 


dx 


(c  ~ c ) 

6 Teq  T 


(46) 


where 


: - f." 

l/2pUZ6  J 


The  equilibrium  value  of  C , (C^  ) was  determined  by  Nash  and  Macdonald, 

eq 

Reference  23,  from  measured  shear  stress  distributions  giving 


eq 


•°25 11 ' if 


(47) 


where  H is  the  local  shape  factor.  Equation  (47)  can  be  expressed  in 


terms  of  the  parameters  u^ , u^  and  U by  evaluation  of  the  integral  rela- 


tions used  to  define  H;  that  is. 


1 
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6 


(48) 


with  the  aid  of  Equation  (43)  . The  equations  for  the  four  unknowns,  u , 

T 

u , <5  and  C , can  be  integrated  once  the  pressure  distribution,  U(x),  and  the 
P T 


streamline  divergence,  1/r (dr/dx),  are  prescribed.  Initial  conditions  are 
obtained  from  the  transition  analysis.  In  this  case,  initial  values  of 
the  momentum  thickness,  0,  and  the  shape  factor,  H,  are  known.  Initial 
values  of  u^,  u^  and  6 can  be  determined  from  the  known  0 and  H using 

Equation  (43).  The  starting  value  for  is  obtained  by  making  the  assump- 


tion that  in  the  region  of  transition, 


C 

T 


The  accuracy  of  the  calculation  method  is  demonstrated  in  Figure  2.  This 
figure  shows  a comparison  between  measured  and  calculated  boundary  layer 
developments  along  a streamline  on  the  U.S.  Airship  Akron.  A further  com- 
parison is  shown  in  Figure  3 for  the  case  of  a boundary  layer  in  a strong, 
adverse  pressure  gradient  approaching  separation.  Particular  reference 
should  be  made  to  the  good  agreement  between  calculated  and  measured  skin 
friction  coefficients. 


ROUGHNESS  AND  AREA  SUCTION 


The  contribution  of  such  roughness  elements  as  rivets  to  the  overall  drag 
of  a helicopter  fuselage  is  small  when  compared  to  drag  due  to  separation. 
However,  as  separation- free  designs  gain  in  interest,  attention  to  such 
items  as  roughness  drag  becomes  important.  The  calculation  method  of  Refe- 
rences 25  and  26  is  used  when  rough  surfaces  are  present  to  determine  the 
effect  of  roughness  on  the  downstream  boundary  layer  development  and  on 
skin  friction  drag.  This  method,  which  is  an  extension  of  Head's  original 
entrainment  method,  is  based  on  the  momentum  integral  equation  (Equation 
(22))  and  the  entrainment  equation  (Reference  27). 


d(<$  - 6 *) dx 


F (H  ) 
1 


(6  - 6*) 


du/dx/u  + dr/dx/r 


(49) 


The  parameter,  H = (6  - 6*)/0,  is  related  to  the  conventional  form  para- 
meter, H,  by  1 

H = G (H)  (50) 

1 
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The  entrainment  functions,  FCH^  and  G (H) , have  been  defined  numerically  as 


and 


G (H) 

for  H < 1.6  or 


F (H  ^ ) = exp  -3.512  - 0.617  An  (Hj  - 3.)  j 

3.3  + exp  0.4667  - 2.722  An  (H  - 0.6798) 


[° 


G (H)  = 3.3  + exp  0,4383  - 3.064  An  (H  - 0.6798) 


] 


(51) 


(52) 


(53) 


for  H > 1.6. 


The  effect  of  roughness  enters  through  the  skin  friction  law  which  can  be 
expressed  as 


(2/c  ) 2 = 5.6  log  uS*/V  + 4.8  - Au  /U  + Au  /U 

f 10  1 T 2 T 


(54) 


where  Au^/U^  represents  the  modifications  of  the  velocity  profile  near  the 

wall  in  the  vicinity  of  roughness  elements.  Au^/U^  has  been  shown  to  be 

a function  of  the  roughness  spacing  (density)  as  well  as  the  roughness 
height.  In  the  density  range,  1 <_  X 5,  where  A is  the  ratio  of  total 
surface  area  to  roughness  area,  the  variation  of  Au^/U.^  with  roughness 

can  be  specified  by 


Au  /U  = 5.o  log  0 k/V  + 17.35  (1.625  log  A -1) 

l l 10  T io 

and  in  the  range  A > 5 by 

Au  /u  = 5.6  log  U k/V  - 5.95  (1.103  log  A - 1) 

1 T 10  T 10 


(55) 


(56) 


The  function  Au^/U^  represents  the  effect  of  pressure  gradient  on  the 
skin  friction,  both  on  smooth  or  rough  surfaces,  and  is  defined  as  follows. 


AU  /U  = 1.253  (G  - 6.7) 

2 l 

for  G ^ 6.7  (adverse  pressure  gradient) , and 


AU  /UT  = 0.404  (G  - 6.7) 

2 i 


(57) 


(58) 
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for  G <6.7  (favorable  pressure  gradient). 


The  parameter,  G,  is  Clauser's  pressure  gradient  parameter  given  by 


G = (2/Cf)  (H  - 1)/H 


(59) 


The  entrainment  equations,  Equations  (49)  and  (50),  may  be  integrated 
simultaneously  with  the  momentum  integral  equation,  with  the  skin  friction 
coefficient  being  determined  at  each  step  in  the  calculation  from  Equation 
(54)  to  give  the  turbulent  boundary  layer  development.  The  calculation  is 
begun  at  the  transition  point  with  initial  values  of  the  shape  factor,  H, 
and  the  momentum  thickness,  0,  determined  from  the  laminar  and  transition 
calculations.  The  roughness  height  distribution,  k/c,  the  roughness  den- 
sity, A,  the  pressure  distribution,  and  the  stream  divergence  distribution 
are  all  assumed  known. 

The  effect  of  area  suction  on  the  development  of  the  turbulent  boundary 
layer  is  readily  accounted  for  by  adding  the  ratio  of  suction  velocity  to 
free-stream  velocity  (Vw/U)  to  the  right-hand  side  of  both  the  momentum 
integral  and  entrainment  equations.  A different  skin  friction  law 
must,  however,  be  used  to  determine  the  effect  of  suction  on  the  local 
skin  friction  coefficient.  This  is  provided  by  the  work  of  Thompson, 
whereby  a skin  friction  law  of  the  form 

Cf  = f(H,  Rq,  V^/U)  (60) 


is  derived.  For  a given  suction  distribution,  the  boundary  layer  devel- 
opment is  calculated  using  Equation  (22)  in  conjunction  with  Equations 
(49)  and  (50) . 

VISCOUS/INVISCID  INTERACTION 


The  effect  of  boundary  layer  displacement  on  the  potential  flow  is  simula- 
ted by  distributing  sources  of  known  strength  on  the  panels  used  to  des- 
cribe the  fuselage  geometry.  The  strengths  of  these  sources  are  deter- 
mined directly  from  the  boundary  layer  solutions  as  q.  = d_  (U.<5.*), 

ds 

where  lb  is  the  streamwise  potential  flow  velocity  at  the  edge  of  the 

boundary  layer,  and  <5^*  is  the  streamwise  displacement  thickness.  The 

addition  of  this  source  distribution  modifies  the  normal  velocity  at  the 
control  point  of  panel  i.  Consequently,  the  boundary  condition  (Refer- 
ence 3)  is  modified  as  follows 

N 

V . = R.  + q.  + V*  a.  . O.  (61) 

ru  i i 13  3 

3 = 1 
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T 


Since  is  known  for  each  iteration,  only  the  right-hand  side  of  Equation 
(61)  is  altered,  giving  (in  matrix  notation) 

[ Aij  ] °j  = -Ri  - qi  <62> 

Because  the  original  geometry  is  not  modified  by  the  use  of  distributed 
sources,  the  aerodynamic  influence  coefficient  matrix  need  not  be  recalcu- 
lated. Subsequent  iterations  between  the  potential  flow  and  boundary  layer 
calculations  result  in  convergent  solutions.  The  alternative  procedure 
(modifying  the  geometry  directly  by  addition  of  the  displacement 
thickness)  while  quite  widely  used  in  two  dimensions,  becomes  untenable  in 
three  dimensions.  Primarily,  this  is  a result  of  having  to  calculate  and 
invert  a new  aerodynamic  influence  coefficient  matrix  at  each  iteration 
due  to  the  change  in  geometry.  Without  considering  the  additional  cost  of 
smoothing  the  geometry  and  redefining  the  panels  at  each  iteration,  the 
additional  cost  of  calculating  and  inverting  the  influence  coefficient 
matrix  at  each  iteration  would  be  prohibitive. 

CALCULATION  PROCEDURE 


The  computer  program  is  made  up  of  a series  of  overlays,  as  shown  in 
Figure  4 . The  executive  program,  DRAG,  controls  the  overall  analysis  by 
calling  in  turn  the  overlays  containing  the  potential  flow  and  boundary 
layer  calculation  methods.  The  calculation  sequence  is  outlined  as 
follows . 

1.  The  input  geometry,  as  represented  by  a series  of  planar  panels, 
is  lofted  in  Program  WBPAN.  An  alternate  slot  is  provided  under 
the  heading  Program  GEOMETRY  for  such  time  as  an  improved  geo- 
metry program  becomes  available. 

2.  The  potential  flow  pressure  field  is  computed  for  the  fuselage 
configuration  in  Program  WBAERO . 

3.  Up  to  six  streamlines  are  calculated  in  Program  STREEM.  These 
streamlines  pass  through  the  midpoints  of  prescribed  panels. 

4.  The  laminar  and  turbulent  boundary  layer  developemnts  are  deter- 
mined for  each  streamline  using  pressure,  arc  length, and  stream- 
line divergence  distributions  calculated  in  Program  STREEM. 
Transition  or  laminar  separation  and  turbulent  separation  are 
predicted,  if  present.  Program  INTGRAL  is  used  for  the  integral 
boundary  layer  analysis;  and  currently,  three  options  are  avail- 
able for  the  turbulent  boundary  layer  analysis;  for  boundary 
layers  on  smooth  surfaces,  rough  surfaces,  or  with  area  suction 


33 


OVERLAY  ( 


Figure  4.  Overlay  Structure  of  the  Calculation  Procedure. 


boundary  layer  control,  A slot  is  provided  for  Program  INSPAN, 
the  finite  difference  method  for  tangential  blowing  boundary 
layer  control.  This  program  will  be  made  available  at  a later 
date. 

5.  Source  distributions  representing  the  boundary  layer  displacement 
effect  are  determined  for  each  streamline. 

6.  The  separated  flow  region  is  determined  in  Program  WAKECON . The 
separation  model  is  prescribed,  and  a new  potential  flow  solution 
(with  the  boundary  layer  source  distribiton  included)  is  obtain- 
ed. Base  pressures  are  determined,  and  the  drag  is  calculated. 

Steps  3 through  6 are  repeated  a second  time  if  required.  The  calculation 
procedure  is  illustrated  in  Figure  5. 
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Figure  5.  'low  Diagram  of  the  Calculation  Procedure. 
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DEVELOPMENT  OF  THE  SEPARATED-FLOW  MODEL 


I 


DESCRIPTION  OF  THE  REAL  FLOW 

Four  regions  are  identifiable  in  the  flow  field  around  a body  with  separa- 
tion (Figure  6)  . 

Region  1 - The  Potential  Flow  region 

The  region  exterior  to  the  boundary  layer  and  separated  wake  is  almost  pre- 
cisely irrotational  since  the  shear  is  so  low  that  viscous  stresses 
impart  a negligible  rotation  to  the  fluid  (originally  irrotational  up- 
stream of  the  body) . All  of  cne  regions  are  almost  solenoidal*  since  we 
are  assuming  the  Mach  number  to  be  low  enough  to  make  compressibility  negli- 
gible. Thus  the  region  is  very  nearly  a potential  flow  (i.e.,  irrotational 
and  solenoidal) . 

Region  2 - The  Boundary  Layer 

The  thin  flow  region  next  to  the  airfoil  surface  has  highdiear,  and  hence, 
viscous  stresses  which  create  significant  vorticity. 

Region  3 - The  Free  Shear  Layer 

The  thin  flow  region  fed  by  the  separating  boundary  layer  has  rotation  but 
only  moderate  shear.  The  vorticity  transport  is  predominantly  by  convec- 
tion, although  diffusion  is  not  insignificant. 

Region  4 - The  Wake 

The  wake  between  the  two  shed  boundary  layers  is  a region  with  low  vorti- 
city and  insignificant  viscous  stresses.  This  region  has  a lower  total 
head  than  that  of  the  onset  flow. 

TWO-DIMENSIONAL  METHOD 

A simple  theoretical  model  of  the  real  flow  is  shown  in  Figure  7.  This 
model  was  developed  in  two-dimensional  flow  before  proceeding  to  the  three- 
dimensional  case.  In  order  to  be  compatible  with  the  WBAERO  method,  the 
body  is  represented  by  a number  of  constant-source  panels,  and  the  free 
shear  layers  by  uniform  vorticity  panels  aligned  in  the  free-stream  direc- 
tion. (Later  developments  should  consider  calculating  the  vortex  panel 
orientation.)  The  vorticity  panels  must  be  attached  to  the  body  at  the 
source  panel  edges  nearest  to  the  separation  point.  This  restriction  is 
necessary  to  avoid  numerical  problems  associated  with  the  abrupt  change 

‘Solenoidal  means  vanishing  divergence,  the  mass  conservation  requirement 
for  a steady,  incompressible  flow. 


REGION  1 


REGION  ] - POTENTIAL  FLOW  REGION 
REGION  2 - BOUNDARY  LAYER 
REGION  3 - FREE  SHEAR  LAYER 
REGION  4 - WAKE 


Figure  6.  Regions  in  the  Real  Flow. 


from  source  to  vortex  panels. 


The  boundary  condition  of  zero  normal  velocity  is  specified  at  a control 
point  in  the  middle  of  each  surface  source  panel.  In  addition,  the  bound- 
ary condition  of  the  streamwise  flow  is  specified  on  the  vortex  panels  at 
points  directly  above  the  source  panel  control  point  downstream  of  the  sepa- 
ration point  (Figure  7 ).  Considering  a symmetrical  case,  we  can  use,  say, 
N source  panels  and  one  vortex  panel  on  half  of  the  body.  This  gives  N + 1 
control  point  equations  in  N unknown  source  panel  densities,  and  one  vorti- 
city  value.  The  equations  are  of  the  form 


Ajk  °k 


V 


C . 
3 


k = 1,  2 N (63) 

j = 1,  2,  . . . , N + 1 


where  A is  the  influence  coefficient  for  the  normal  component  of  veloc- 
3 ^ 

ity  induced  at  the  j^1  panel's  control  point  by  the  k source  panel;  Eh 

is  the  influence  coefficient  for  the  normal  component  of  velocity  induced 


at  the  panel's  control  point  by  the  vortex  panel;  and  C.  is  the  normal 

component  of  the  onset  flow  at  the  panel’s  control  point.  Later,  when 

the  boundary  layer  procedure  is  included  to  calculate  the  separation  point, 
CL  will  include  the  boundary  layer  source  distribution  representing  dis- 
placement effect  (see  Section  on  Viscous/Potential  Flow  Interaction)  • 


Solution  of  the  system  of  linear  equations  gives  the  source  distribution 
and  vorticity  value  from  which  the  surface  pressures  are  evaluated.  In 
the  attached  flow  region,  the  pressure  coefficient  is 

2 


In  the  separated  flow  region,  there  is  a reduction  in  total  head  from  that 
of  the  free  stream.  Here,  therefore,  the  pressure  coefficient  is 


The  reduction  in  total  head,  AH,  can  be  evaluated  by  considering  the  condi- 
tions on  the  free  shear  layer.  On  the  outer  edge  of  the  layer 
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and  on  the  inner  edge 


V. 

x 


2 


c 

pi 


1 


Ah 


2 


and  since  the  layer  cannot  support  a pressure  difference  across  it,  then 
C = C , and 

P„  P,- 


or  AH  = Jjp(V  - V.)  (V  + V.) 

r o 1 o 1 


But  (V  - Vk)  is  the  vorticity  value  in  the  layer,  i.e.,  y,  which  is  part 
of  the  solution,  and  ~ Vj  is  the  mean  velocity,  V^,  say,  in  the  layer. 

Therefore , 


AH  = pVMy 


and  the  pressures  in  the  separated  region  are 


2V 


The  mean  velocity,  V , in  the  shear  layer  is  readily  calculated  once  the 
singularity  strengths  are  known. 

Initial  computations  based  on  a circular  cylinder  showed  that  the  peak 
suction  upstream  of  separation  was  underestimated  when  using  a long  wake 
panel.  It  was  found  that  a shorter  wake  length  gave  good  agreement  with 
experimental  pressures  over  the  whole  cylinder.  This  sensitivity  of  the 
solution  to  the  assumed  wake  length  is  illustrated  in  Figure  8 for  a 
circular  cylinder  with  separation  at  10511.  The  wake  length  which  gives 
good  agreement  with  the  experimental  base  and  peak  suction  pressures  (from 
Reference  28)  in  this  case  is  about  20%  larger  than  the  cylinder  diameter. 
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Figure  8.  Base  Pressure  and  Peak  Suction  on  a Circular  Cylinder  as  a 
Function  of  Wake  Panel  Length. 


Since  the  "correct"  wake  length  is  not  generally  known  a priori,  the  problem 
is  now  a nonlinear  one.  however,  it  was  observed  that  the  "correct"  wake 
length  produced  essentially  a constant  pressure  distribution  along  the  outer 
edge  of  the  wake  panel.  It  was  possible  to  make  the  problem  linear  again  by 
using  the  constant  pressure  criterion  combined  with  a multiple  segment  vor- 
tex panel.  This  required  additional  control  points  and  equations.  Although 
satisfactory  solutions  were  obtained  directly,  i.e.,  without  iteration,  it 
was  felt  that  the  procedure  might  not  be  general  enough  for  application  to 
other  shapes  since  some  wake  length  was  still  prescribed.  An  iterative 
technique  was,  therefore,  developed  and  tested  on  the  circular  cylinder 
case.  This  technique  used  the  criterion  of  zero  tangential  velocity  dif- 
ference between  two  points  on  the  wake  panel  (Figure  9)  . Over  a short  dis- 
tance along  the  wake,  the  tangential  velocity  increases  if  the  wake  is  too 
short,  and  it  decreases  if  the  wake  is  too  long.  Evaluating  the  tangen- 
tial velocity  increment  for  two  wake  length  assumptions  (e.g.,  1.  and  1.5 
times  the  base  diameter)  enables  a wake  length  to  be  interpolated,  which 
gives  zero  velocity  increment.  The  process  can  be  continued,  but  three 
iterations  are  generally  sufficient  (Figure  10) . 

The  "correct"  wake  length  also  gives  a small  plateau  in  the  normal  velocity 
component  distribution  (Figure  9 ) . This  implies  that  the  choice  of  the 
control  point  location  on  the  wake  panel  is  not  critical. 

Figure  11  shows  the  calculated  pressure  distribution  on  the  circular  cylin- 
der, using  two  panel  densities,  compared  with  experimental  results  for  a 
Reynolds  Number  of  8.4  x 106  (Reference  28).  The  separation  angle,  i.e., 
105°,  was  prescribed  at  this  stage  in  the  model  development.  The  wake 
length  iteration  was  stopped  after  the  third  pass . The  20-panel  case  (on 
half  the  cylinder)  gave  very  good  agreement  with  the  experimental  pressure 
distribution,  particularly  in  the  separated  region,  but  there  is  an  "over- 
shoot" in  pressure  just  upstream  of  the  separation  point. 

The  case  with  only  10  panels  also  gave  excellent  agreement  with  experiment 
in  the  separated  flow  region  pressures,  but  it  underestimated  the  suctions 
upstream  of  the  separation  point. 


THREE-DIMENSIONAL  METHOD 


Testing  the  Potential  Flow  Model 

The  two-dimensional  procedure  described  in  the  previous  subsection  was  ex- 
tended to  the  three-dimensional  form  and  applied  to  a sphere.  In  this 
case,  the  wake  vorticity  panels  form  a cylinder  behind  the  sphere.  Initial 
calculations  showed  a different  behavior  in  the  iterative  wake-length  pro- 
cedure from  that  shown  in  the  two-dimensional  case.  Therefore,  the  effect 
of  wake  length  on  base  pressure  was  reevaluated,  and  the  results  are 
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shown  in  Figure  12  for  three  prescribed  separation  angles.  The  base  pres- 
sure now  appears  to  have  an  asymptotic  value  as  wake  length  increases.  It 
is  virtually  insensitive  to  wake  length  for  lengths  above  about  twice  the 
separated  wake  diameter.  The  iterative  wake  length  procedure,  which  was 
found  to  be  necessary  in  the  two-dimensional  case,  is  evidently  not  re- 
quired in  the  three-dimensional  case. 


Figure  13  shows  the  variation  of  the  calculated  base  pressure  on  the  sphere 
with  a prescribed  separation  angle.  The  calculated  values  are  enclosed 
lines  representing  wake  lengths  of  one  separated-wake  diameter  and  essen- 
tially infinity.  The  enclosed  region  is  not  very  wide,  particularly  at  the 
larger  separation  angles  (i.e.,  higher  Reynolds  numbers) . Two  experimental 
points  are  included  from  Reference  29  for  Reynolds  numbers  .425  x 106  and 
.157  x 10' . In  view  of  the  difficulty  in  obtaining  the  separation  point 
from  the  experimental  results,  the  agreement  between  the  theoretical  model 
and  experiment  is  remarkably  good. 

The  rapid  variation  in  base  pressure  with  separation  angle.  Figure  13,  will 
pose  a problem  with  the  present  theoretical  model  in  that  the  necessary 
displacement  of  the  wake  panels  from  the  separation  line  to  the  nearest 
source  panel  edge  will  cause  a change  in  predicted  base  pressure.  This 
might  eventually  lead  to  a requirement  of  a more  flexible  model  that  is 
not  constrained  to  follow  panel  edges. 

FEATURES  OF  THE  COMPLETE  PROCEDURE 

The  full  potential  flow  . iscous  flow  iterative  procedure  described  earlier 
was  investigated  next.  The  separation  points,  instead  of  being  prescribed 
as  in  the  previous  subsection,  are  now  caiculated  in  the  boundary  layer 
procedure  based  on  the  previous  potential  flow  pressure  distribution  along 
calculated  streamlines.  Preliminary  calculations  highlighted  several  prob- 
lems which  required  scan  i ti  ns  t b<  applied  to  the  method.  These 

are  described  below. 


"Rectangular"  Paneling  Arrangement 

The  paneling  of  the  surface  containing  the  separated  region  must  form  a rec- 
tangular grid  at  this  stage;  i.e.,  the  number  of  rows  in  each  column  must 
be  constant  from  the  start  of  the  body  to  a station  well  downstream  of  the 
expected  separation  line.  This  restriction  is  necessary  at  this  stage  for 
the  streamline  calculation  procedure  and  for  the  procedures  which  find  the 
"separation"  panels  and  which  redistribute  the  boundary  layer  sources  and 
skin  friction  drag.  The  procedure  that  subtracts  Ah  from  pressures  on 
panels  downstream  of  the  separation  line  in  the  "rectangular"  paneling 
region  also  relies  on  this  restriction.  It  would  be  possible  to  develop 
more  general  procedures  to  take  advantage  of  the  arbitrary  paneling  faci- 
lity of  the  WBAERO  program,  but  such  procedures  would  be  considerably  more 
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SPHERE  DIAM.  3 1 


Figure  12.  Effect  of  Wake  Length  on  Calculated  Base  Pressure  on 
a Sphere  for  Three  Separation  Angles. 


complicated.  At  this  stage  of  the  development  of  the  wake  model,  the  simple 
"rectangular"  scheme  is  adequate. 

The  Separation  "Line" 

The  wake  paneling  program  finds  the  surface  panels  that  have  their  upstream 
edge  nearest  to  the  predicted  separation  points.  It  then  finds  the  columi 
of  panels  that  has  the  most  separation  panels  defined  in  it  and  uses  the 
upstream  edge  of  that  column  as  the  starting  point  for  the  wake  vorticity 
panels.  Each  surface  source  panel  in  that  column  has  a wake  vorticity  panel 
iched  to  its  upstream  edge.  This  restriction  was  found  to  be  necessary 
because  breaks  in  the  vorticity  paneling  give  distorted  vorticity  solu 
tions  and  poor  convergence  characteristics  in  the  iterative  sol  ;tion  pro- 
cedure . Such  breaks  occur  when  the  vorticity  paneling  is  allowed  to  "jump  " 
from  one  column  to  another  in  order  to  represent  a separation  line  that  is 
inclined  relative  to  the  paneling.  Tnis  restriction  can  result  in  a i >or 
representation  of  the  separation  line,  particularly  in  angle  of  attack 
cases.  The  problem  would  only  be  removed  by  allowing  the  separation 
panels  to  be  attached  anywhere  on  tin  surface  panels.  Such  a feature  would 
not  be  possible  with  the  present  constant  source/vortex  panel  model  because 
of  the  interaction  between  the  panel  edges  and  the  control  points  and  be- 
cause of  the  abrupt  change  in  the  type  of  singularity. 

Boundary  Condition  on  the  Wak»  _ J\.  ■ ■ h 

The  vorticity  level  on  the  wake  panels  is  determined  after  th<  boundary 
dition  of  streamwi.se  flow  is  added  on  the  .-ontr-.l  points  near  the  up:  t si 
edge  of  the  panels.  This  is  analogous  to  a Kutta  trailing-edge  condition, 
and  it  can  force  the  local  flow  to  deviate  appreciably  from  the  real  flow 
direction.  This  is  a possible  contributing  factor  to  the  low  calculated 
peak  suction  just  upstream  of  th<  separation. 

Solution  of  the  Singularity  Strengths 

Because  of  the  large  number  of  panels  needed  to  represent  a body,  iterative 
techniques  are  favored  for  the  solution  of  the  source  i;  1 vortex  strengths. 
Unfortunately,  because  of  adverse  interference  between  source  and  vortex 
panels  near  the  separation  line,  iterative  techniques  failed  to  converge. 

A modified  technique  was  used,  therefore,  which  ombines  iterative  and 
direct  solutions.  This  is  basically  the  same  as  the  procedure  described 
in  Reference  3;  but  the  direct  solution  part  covers  the  vorticity  and 
part  of  the  source  partitions,  so  that  the  main  source  vortex  interference 
terms  are  included.  The  procedure  starts  with  the  Jauss-Seidel  method 
with  the  initial  assumption  of  zero  source  strength  and  0.1  for  the  vorti- 
city. (The  vorticity  value  is  roughly  equal  to  the  free-stream  velocity 
since  it  represents  the  jump  in  tangential  velocity  from  approximately 
free  stream  to  approximately  zero.  The  solution  evaluates  y/4ir.)  At  a 
certain  point  in  the  matrix,  the  Gauss-Seidel  procedure  is  stopped  and 
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a direct  solution  performed  for  the  remaining  equations.  The  basic  equa- 
tions are  of  the  form: 


T 


N 


E A^k 


M 


E v 


N + M 


This  differs  from  the  two-dimensional  form  (Equation  (63) ) in  that  there  are 
now  M wake  panels  instead  of  1.  If  the  Gauss-Seidel  operation  is  taken  up 
to  j = m,  say,  then  for  the  rest  of  the  equations  we  can  write 

EV»  ‘ EsJk°k+  * S; 

k=l  k=m+l  k=l 

j =m+  1,  m+  2 , ...»  N+M 


Using  the  present  Gauss-Seidel  solution  for  0 , k = 1,  2,  ...»  m,  the 

first  term  can  be  evaluated  and  combined  with  the  right-hand  side  for  all 
of  the  remaining  equations.  This  leaves  a square  system  of  equations  which 
is  solved  by  a direct  method.  The  direct  solution  for  o , k = i + 1,  m + 2, 

. ..,  N,  and  for  k = 1,  2,  ...»  M is  then  used  for  the  next  pass  through 

the  Gauss-Seidel  procedure  and  so  on.  Currently,  the  direct  solution  is 
performed  for  the  last  140  equations,  which  adequately  covers  the  region  of 
close  interaction  between  the  source  and  vortex  panels.  The  procedure 
typically  converges  within  six  iterations  for  a residual  limit  of  10  3 . If 
the  number  of  unknowns  is  less  than  141,  then  the  whole  solution  is  by  the 
direct  method.  The  attached  flow  solution,  which  involves  source  panels  only, 
is  always  obtained  by  the  unmodified  Gauss-Seidel  technique. 

Boundary  Layer  Source  and  Skin  Friction  Distr ibut ion 

Boundary  layer  source  values  and  skin  friction  values  are  known  only  along 
each  streamline.  Because  the  panel  history  is  known  along  each  streamline, 
the  source  and  skin  friction  values  can  be  set  on  those  panels  crossed  by 
streamlines.  The  values  on  other  panels  are  set  by  interpolation  down  each 
column.  Downstream  of  the  separation  line,  the  boundary  layer  source  and 
skin  friction  values  are  set  to  zero. 
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CALCULATION  AND  DISCUSSION  OF  RESULTS 


SPHERE  CASE 

A sphere  calculation  using  112  panels  was  performed  for  a Reynolds  Number 
of  .425  x 106.  The  resulting  pressure  distribution  for  two  iterations  is 
shown  in  Figure  14  together  with  an  experimental  distribution  from  Refer- 
ence 29.  Also  included  are  the  exact  and  calculated  attached  flow  pres- 
sure distributions.  The  first  iteration  result  shows  a later  separation 
angle  than  the  experiment  and  a higher  base  pressure;  i.e.,  C is  .3 

PBASE 

instead  of  .1.  The  second  iteration  shows  an  earlier  separation  than  the 
experiment  and  a lower  base  pressure;  i.e.,  is  .02  instead  of  .1.  Fur- 
ther iterations  were  not  attempted  at  this  stage  because  it  is  suspected 
that  earlier  calculated  separation  is  being  forced  after  the  first  itera- 
tion because  of  the  positive  pressure  overshoot  just  upstream  of  the  wake 
model  vorticity  panels  (Figure  14) . This  augments  the  adverse  pressure 
gradient  over  the  rear  of  the  body.  Both  iterations  show  a loss  in  peak 
suction  on  the  sphere  (Figure  14) . A factor  which  may  be  causing  the  re- 
duced peak  suction  and  also  the  positive  pressure  overshoot  is  the  boundary 
condition  that  forces  the  flow  to  be  stream-wise  just  after  separation. 

An  improved  boundary  condition  for  the  wake  panels  might,  therefore,  be 
necessary . 


BO 10 5 FUSELAGE  CASES 

The  method  was  applied  to  the  B0105  fuselage  using  two  paneling  arrange- 
ments. The  first  case  had  72  panels.  The  "rectangular"  paneling  region 
had  6 rows  by  10  columns  and  covered  the  main  part  of  the  fuselage  and  the 
start  of  the  boom  (Figure  15) . The  second  case  had  241  panels  (Figure  16) 
with  9 rows  by  25  columns  on  the  "rectangular"  paneling  area.  The  calcula- 
ted pressure  distributions  along  waterline  6 (Figure  16)  for  the  two  panel- 
ing schemes  are  compared  in  Figure  17  for  the  second  potential  flow/viscous 
flow  iteration.  This  clearly  demonstrates  the  need  for  fairly  dense  panel- 
ing to  define  the  peak  suction  areas  adequately.  Both  distributions  show 
the  positive  pressure  overshoot  just  upstream  of  the  separation  that  was 
observed  in  the  sphere  and  cylinder  cases. 

Figure  18  shows  the  calculated  streamlines  from  the  dense  paneling  case 
for  two  iterations.  The  starting  points  for  the  streamline  calculations 
were  in  the  column  of  panels  at  x-station  36.5.  The  changes  in  the  stream- 
line paths  are  fairly  small,  and  the  calculated  separation  points  on  each 
streamline  are  in  good  visual  agreement  with  experimental  observations 
(Reference  2) . 
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STREAMLINES  CALCULATED  SEPARATION  POINTS 
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Figure  18.  Calculated  Streamlines  and  Separation  Points  on  the  B0105  Fuselage  for 
Iterations  (241-Panel  Case). 


The  starting  points  for  the  wake  vortex  panels  (Figure  18)  show  a reason- 
able representation  of  the  separation  line  for  the  first  iteration,  except 
on  the  top  surface  of  the  body  where  the  calculated  separation  line  moves 
downstream.  However,  on  the  second  iteration,  there  is  an  appreciable  for- 
ward displacement  of  the  wake  paneling  from  the  calculated  separation  line, 
particularly  on  the  upper  part  of  the  body.  This  would  imply  a shift  in 
the  predicted  base  pressure,  following  the  results  observed  in  Figure  13. 

Calculated  and  experimental  (Reference  2)  pressure  distributions  are  com- 
pared along  three  lines:  top  centerline,  waterline  6 (see  Figure  18),  and 

bottom  centerline,  in  Figure  19  (a) , (b)  and  (c) , respectively.  The  calcu- 
lations are  from  the  241-panel  case  and  the  attached  potential  flow,  the 
first  and  second  potential  flow/viscous  flow  iteration  results  are  presented. 
The  first  iteration  did  not  predict  separation  on  the  top  of  the  body  within 
the  "rectangular"  paneling  area,  but  the  separation  paneling  from  the  lower 
part  of  the  body  was  continued  onto  this  region  to  avoid  an  abrupt  change 
in  the  vortex/source  paneling.  On  the  second  iteration,  separation  was 
predicted  at  station  42.8  on  the  top  streamline;  but,  to  avoid  breaks  in 
the  wake  model,  the  vortex  paneling  started  at  station  38  to  suit  the  major- 
ity of  the  predicted  separation  points.  The  effect  of  this  forward  "sepa- 
ration" is  reflected  in  the  top  centerline  pressure  distribution  in  Figure 
19(a).  The  good  agreement  with  the  experimental  base  pressure  on  this  line 
is,  perhaps,  fortuitous  in  view  of  the  displaced  wake  model.  As  in  the 
sphere  case  (figure  14) , there  is  a positive  pressure  overshoot  and  an 
underprediction  of  peak  suction  just  ahead  of  the  separated  region.  Up- 
stream of  this  region,  the  pressures  are  in  good  visual  agreement  within 
the  experiment. 

The  waterline  6 pressure  distribution,  Figure  19(b)  , shows  similar  tenden- 
cies to  the  top  centerline  distribution,  but  in  this  case  the  base  pressure 
on  the  second  iteration  is  about  .17  too  high.  This  is  not  consistent  with 
the  fact  that  the  wake  panels  are  forward  of  the  predicted  separation  line 
(see  Figure  13).  Figure  19(c)  shows  the  bottom  centerline  distributions. 
Here,  the  base  pressure  at  the  second  iteration  is  in  good  agreement  with 
the  experiment.  The  separation  paneling  in  this  case  is  closest  to  the 
calculated  separation  points  and  the  pressure  distributions  compare  reason- 
ably well  except,  as  observed  before,  the  calculations  underestimate  the 
rear  peak  suction. 

The  calculated  positive  pressure  overshoot  just  upstream  of  the  separation 
paneling  is  a detrimental  feature  in  all  three  distributions.  Not  only 
does  this  automatically  cause  earlier  separation  to  be  predicted  on  the 
next  iteration  (because  of  the  more  adverse  pressure  gradient) , but  it  has 
a large  influence  on  the  integrated  pressure  drag.  The  calculated  drag 
breakdown,  D/q,  at  the  second  iteration  is  as  follows  (in  square  feet) : 
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Figure  19.  Continued. 


BO  105  TOP  CENTERLINE 


(c)  Bottom  Centerline 


Skin  Friction  .045 

Integrated  Pressure  -.087 

Total  -.042 


The  positive  pressure  overshoot  is  estimated  to  have  lost  approximately 

rated'TL  .°f  °/q-  Base  Pressure  and  peak  suction  errors  are  esti- 
mated to  have  lost  approximately  .07  and  .04  square  foot  respectively. 

The  boom  gives  additional  small  drag  errors  because  it  is  not  fully  Leafed 

S‘«»-  f"  subtracted  fro.  all  panel,  L.Hf 

the  wake  (x.e.,  not  just  those  in  the  "rectangular"  paneling  area),  but 
this  would  require  an  extensive  search  procedure  to  check  the  geometry  of 
each  panel.  1 
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CONCLUSIONS 


1.  The  separated  flow  model  has  demonstrated  reasonable  accuracy  in  the 
predicted  surface  pressures  upstream  of  the  separation  line  with  fair 
agreement  in  base  pressure.  However,  improvements  are  needed  in  the  rear 
peak  suction  regions  and  the  positive  pressure  overshoot  just  upstream 
of  the  separation  line  should  be  removed. 

2.  Certain  improvements  are  needed  in  the  model  to  allow  more  flexibility 
in  the  wake  paneling.  A more  accurate  representation  of  the  separa- 
tion line  is  required,  rather  than  moving  the  wake  panels  to  the  nearest 
source  panel  edges . Such  an  improvement  should  result  in  more  accurate 
predictions  of  base  pressure. 

3.  The  calculated  streamlines  and  separation  points  are  in  good  visual 
agreement  with  experimental  observations. 
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SUGGESTED  PROGRAM  IMPROVEMENTS 


The  detrimental  features  of  the  source/vortex  panel  model  can  be  largely 
removed  by  going  to  a vorticity  panel  method.  Such  a method  has  been 
developed  very  recently  for  predicting  lift  coefficients  beyond  the  stall 
in  two-dimensional  flow.  Applied  to  the  drag  prediction  problem,  the  vorti- 
city model  would  have  surface  panels  of  linearly  varying  vorticity  (Figure 
20).  The  wake  vorticity  panels  would  be  similar  to  those  in  the  present 
work,  but  the  vorticity  level  would  be  the  same  as  the  value  on  the  surface 
panel  at  the  attachment  point.  The  "Kutta"  condition  can  be  applied  which 
makes  the  lower  vorticity  panel  to  be  equal  but  opposite  to  the  upper  one; 
this  feature  would  extend  the  present  method  to  allow  lifting  cases  to  be 
considered.  Inside  the  separated  wake  region,  vorticity  is  constrained  to 
be  zero  just  downstream  of  the  separation  points;  the  rest  of  the  corner 
vorticity  values  in  between  are  evaluated  in  the  solution.  Boundary  con- 
ditions, therefore,  are  applied  only  on  the  body  surface;  hence,  the  wake 
panel  boundary  points  of  the  source /vortex  model  are  eliminated.  This 
feature  contributes  to  an  improved  prediction  of  upstream  suction;  e.g.,  in 
the  cylinder  case,  12  vorticity  panels  give  very  good  comparison  with  the 
experimental  pressure  distribution  (Figure  21) . They  are  equivalent  to  20 
source  panels  and  yet  do  not  show  the  positive  pressure  overshoot.  Positive 
pressure  "overshoots"  and  "undershoots"  have  been  obtained  from  the  vor- 
ticity model  when  using  a large  number  of  panels  on  an  airfoil.  Correct 
wake  orientation  seems  to  be  the  answer  to  this  problem,  and  a relaxed  wake 
procedure  might  be  called  for  eventually. 

With  the  surface  vorticity  model,  the  vorticity  is  continuous,  moving  from 
the  surface  onto  the  wake  panel.  There  are  no  abrupt  changes  in  singula- 
rity strength  or  type.  Arbitrary  wake  paneling  is,  therefore,  possible; 
i.e.,  the  wake  panels  can  oe  attached  along  any  line  on  the  surface  pane’s. 
Accurate  representation  of  the  separation  line  would,  therefore,  be  possible, 
with  this  model.  This  would  remove  the  problems  associated  with  the 
source/vortex  wake  paneling  restrictions.  In  addition,  the  vorticity  con- 
tinuity feature  would  allow  the  separation  modeling  to  finish  in  the  middle 
of  a surface  - the  vorticity  would  simply  stay  on  the  surface  where  sepa- 
ration is  not  predicted.  Because  of  the  continuity,  numerical  solution 
of  the  equations  should  be  better  behaved  than  in  the  case  of  the  source/ 
vortex  model. 
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LINEAR  VORTICITY  PANELS 


CONTROL  POINTS 

CONSTANT- VORTICITY 


ACROSS  SEPARATION  POINT, 
i.e.,  FROM  THE  UPSTREAM 
PANEL  ONTO  THE  WAKE 
PANEL. 


Figure  20.  Separated  Flow  Model  Based  on  Vorticitv  Panels. 
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APPENDIX  I.  PROGRAM  USER’S  GUIDE 


PROGRAM  INPUT  DESCRIPTION 


Card  Set 


Input 


Format 


1. 


Case  Title 


8A10 


2. 

RNB,  TRIPUP,  OPTION,  TRMAX,  XROUGH,  REFC, 
UIN,  XGEM 

8F10.0 

3. 

UKU,  DEN,  XO,  EPS,  XN ROUGH 

7F10.0 

4. 

VRD(l)  ...  VRD (NROUGH) 

7F10.0 

5. 

XPRINT , XSKIP , REFX,  REFZ,  CREF,  PRINT,  CASE 

7F10.0 

6. 

NPANEL,  NBX,  MBX,  NWX,  MWX,  NVX,  MVX,  ILX,  NPTS 

915 

7. 

TEXT  Card:-  GEOMETRY  INPUT 

8A10 

8. 

CASE,  PLOT,  SIM,  ISAVE,  PRINT 

5110 

9A. 

SINGPA,  NOPAN 

2110 

9B. 

X(I)  , Y (I)  , Z(I) 

7F10.0 

9C. 

NPAN  (1)  ...  NP AN (NOPAN) 

7110 

10A. 

NB 

110 

10B. 

XBE,  YBE , ZBE,  MB,  OPT,  FLAG 

3F10.0,  3110 

IOC. 

B ( J)  , A (J)  , D (J) 

7F10.0 

10D . 

D(J) , B ( J) , A (J) 

3F12.0 

11A. 

NW,  KOORD 

2110 

11B. 

XBE,  YBE,  ZBE,  CHRD,  ALF,  XAL,  MW,  OPT,  FLAG, 
DEL 

6F10,  415 

11C. 

DELTA,  YO,  ZO 

7F10.0 

11D. 

B ( J)  , A (J) 

7F10.0 

HE. 

B (J) , D ( J ) , A (J) , C(J) 

5F12.0 

12. 

WAKE,  POINT 

7F10.0 

13A. 

XLP,  YLP , ZLP 

7F10.0 

13B. 

XLP,  YLP,  ZLP 

7F10.0 

14A. 

NV 

110 

14B. 

XBE,  YBE,  ZBE,  MV,  OPT,  FLAG,  SIMPOT 

3F10.0,  4110 

14C. 

B (J) , A ( J ) , C(J) , D (J) 

7F10.0 

14D1 . 

A,  B,  D 

7F10.0 
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26.  Dvorak,  F.A.,  "The  Calculation  of  Compressible  Turbulent  Boundary  Layers 
with  Roughness  and  Heat  Transfer,  AIAA  Journal,  Vol.  10,  No.  11, 

November  1972. 

27.  Head,  M.R. , "Entrainment  in  the  Turbulent  Boundary  Layer",  R and  M 3152, 
1958,  Aeronautical  Research  Council. 

28.  Roshko,  A. , "Analytical  and  Experimental  Studies  on  Drag  and  Flow  of 
Two-Dimensional  Bodies",  NACA  Technical  Notes  2913  and  3169,  1953  and 
1954. 

29.  Fage , A.,  "Experiments  on  a Sphere  at  Critical  Reynolds  Numbers",  Aero. 
Res.  Comm.  Tech.  Rept.  R.&  M.  17566,  September  1936. 


Card  Set 

Input 

Format 

14D2. 

A,  B,  D 

7F10.0 

14D3. 

A,  B,  D 

7F10.0 

15. 

TEXT  Card:-  AERODYNAMIC  CALCULATIONS 

8A10 

16. 

NIT,  I EPS,  I TYPE 

3110 

17. 

COMPT,  SECT 

2110 

18. 

REFA,  REFL,  XOO,  X25 

7F10.0 

19. 

KUT,  NBV,  NV (1 ) ...  NV (5) 

7110 

20. 

KOMPR,  POINTS 

2110 

21. 

NMA 

110 

22. 

MA 

7F10.0 

23. 

NAL 

110 

24. 

ALPHA,  BETA 

7F10.0 

25. 

IS 

110 

26. 

IA,  IE 

2110 

27. 

IREI 

110 

28. 

11(1) , 11(2)  ...  ETC. 

1415 

29. 

DELY,  REFL,  XLE 

7F10.0 

30. 

TEXT  CARD:-  STREAMLINE  CALCULATIONS 

8A10 

31. 

ENDFILE 

110 

32. 

NLIN,  NSP(l),  ...  NSP(NLIN) 

1415 

33. 

6/7/8/9 

71 


Description  of  Input  Variables 


Card  1 - General  Identification  - Card  1 contains  any  dasired  identifyin  ; 
information  in  Columns  1-80. 

Card  2 - General  Parameters 


Column 


1-10 


31-40 


Variable 


RNB 


11-20  TRIPUP 


21-30  OPTION 


41-50  XROUGH 


Value  Description 

arbitrary  Reynolds  number  basec  on  reference  chord 
and  free-streum  velocity 
U^C/V  x (10~6 ) 

N.B.  If  the  input  geometry  is  in  dimen- 
sional units,  the  Reynolds  number  will 
be  input  per  unit,  i.e..  Re/inch. 

" Trip  location  (x/c) 

1.  No  tripping  desired 

0.  Boundary  layer  will  trip  where  specified 

by  TRIPUP 

1.  Program  will  test  on  boundary  layer 
momentum  thickness  at  trip  location. 

If  Rq  < 200,  trip  location  will  be  re- 
positioned to  point  where  Ry  200. 

This  deters  the  user  from  specifying  an 
unrealistically  early  trip  location 

Maximum  number  of  iterations  between  in- 
viscid  and  boundary  layer  modules 

-1.  Standard  boundary  layer  calculation  for 

smooth  surfaces 


TRMAX  arbitrary 


Boundary  layer  calculation  with  area 
suction  (needs  cards  3 and  4) 


51-60 


REFC  arbitrary 


Boundary  layer  calculation  with  surface 
roughness  (needs  cards  3 and  4) 

Reference  chord  in  inches  for  determina- 
tion of  surface  distance  in  INSPAN 


r 


Column 

Variable 

Value 

Description 

61-70 

UIN 

arbitrary 

Eree-stream  velocity  in  feet  per  second 
for  use  in  INSPAN 

71-80 

XGEM 

0. 

Standard  case;  program  will  execute  for 
complete  case 

1. 

( roqram  will  assemble  and  print  out  geo- 
metry variables  only  and  stop 

N.B.  : 

REFC,  UIN  and 

XGEM  currently  not  available. 

Card  3 

- Rouqhness  Parameters 

Column 

Variable 

Value 

Description 

1-10 

UKU 

30.-70. 

Value  of  roughness  Reynolds  number  where 
,'ffect  of  roughness  is  independent  of 
viscosity  (see  Ref.  1) 

11-20 

DEN 

2.-100. 

Density  of  roughness  element  spacing. 
(For  sand  grain  roughness,  use  3.175) 

21-30 

XO 

.005-. 02 

Initial  value  of  skin  friction  coeffi- 
cient (.007  is  usual  starting  value) 

31-40 

EPS 

.0001 

Error  bound  in  skin  friction  calculation 

41-50 

XNROUGH 

arbitrary 

Number  of  input  values  of  suction  velo- 
city, v /U,  or  roughness  heights,  k/c 
w 

Card  4 

- Roughness  or 

Suction  Distribution 

Column 

Variable 

Value 

Description 

1-10 

• 

• 

VRD(l) 

• 

• 

arbitrary 

(floating 

point) 

Roughness  height  or  area  suction  velocity 

(k/c  or  V/U  ) 

00 

61-70 

VRD(N ROUGH) 

«« 

tf 
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Card  5 - General  Parameters 


Column 

Variable 

Value 

Description 

1-10 

XPRINT 

0. 

Suppresses  printing  of  cross-flow  inte- 
gral thickness  from  IBL 

1. 

Extra  printing 

11-20 

XSKIP 

1. 

Print  option  in  INTGRAL,  every  integra- 
tion step  is  printed 

10. 

Every  tenth  step  is  printed 

21-30 

REFX 

arbitrary 

Reference  (x/c)  location  for  calculation 
of  moment  coefficient 

31-40 

REFZ 

It 

Reference  (z/c)  location  for  calculation 
of  moment  coefficient 

41-50 

CREF 

1. 

Redundant 

51-60 

PRINT 

0. 

II 

61-70 

CASE 

1. 

Number  of  angle-of-attack  or  Mach  number 
variations  for  a given  geometry.  Cur- 
rently limited  to  1 

Card  6 - Variable  Dimension  Options 


Column 

Variable 

Value 

Description 

1-5 

NPANEL 

.LE.1500 

Number  of  source  panels  and  vortex 
lattices  on  one  side  of  plane  of  syircnetry 

6-10 

NBX 

2-70 

Number  of  body 

sections 

11-15 

MBX 

3-60 

Maximum  number 
body  section 

of  input  points  for 

any 

16-20 

NWX 

2-40 

Number  of  wing 

sections 

21-25 

MWX 

3-59 

Maximum  number 
section 

of  ordinates  at  any 

wing 
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Column 

Variable 

Value 

Description 

26-30 

NVX 

0-40 

Number  of  streamwise  vortices  in 
lattice 

body 

31-35 

MVX 

2-60 

Number  of  bound  vortices  in  body 

lattice 

36-40 

I LX 

2-35 

Sum  of  wing  and  body  lattices 

41-45 

NPTS 

.LE.1500 

Number  of  off-body  points 

Geometry  Input  Cards  - If  the  configuration  is  symmetrical  about  the  x,z 
plane,  geometrical  input  is  required  for  only  one  side  of  the  configuration. 
The  convention  used  herein  is  to  present  that  half  of  the  configuration 
lying  on  the  positive  y side  of  the  x,z  plane.  If  the  configuration  is  not 
symmetric,  complete  geometrical  input  is  required. 


Card  7 - General  Identification  - Card  7 contains  any  desired  identifying 
information  in  Columns  1-80. 


Column  Variable  Value 


10  CASE  1 

2 

3 

20  PLOT  0 

1 

30  SIM  0 

1 


Description 

Isolated  body  only 

Isolated  wing  only 

Wing-body  combinations 

No  plot  output  (currently  this  is  the 
only  option  available) 

Plot  output  requested 

Conf iguration  symmetric  about  x,z  plane; 
panel  geometry  required  on  one  side  only 
(normal  case) 

Configuration  symmetric  about  x,z  plane. 
Panel  geometry  input  required  on  one 
side  only;  panel  geometry  output  calcu- 
lated for  both  sides.  (Used  when  analy- 
zing symmetric  conf iguration  in  yaw.) 

Not  currently  available 

Asymmetric  configuration.  Panel  geome- 
try input  required  for  both  sides. 
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Column 

Variable 

Value 

Description 

40 

ISAVE 

0 

Geometry  and  influence  coefficient  matri- 
ces not  saved 

1 

Geometry  and  influence  coefficient  matri- 
ces saved  in  previous  run  to  be  used 
(TAPE  11  must  be  requested) . This  option 
not  yet  checked  out 

-1 

Geometry  and  influence  coefficient  matri- 
ces to  be  saved  on  TAPE  11 

50 

PRINT 

0 

Normal  output 

1 

Optional  output  1 . Includes  panel  geo- 
metry, coordinate  transformation  matri- 
ces, and  panel  forces  and  moments 

2 

Optional  output  2.  Panel  velocity  com- 
ponents and  influence  coef f icients . 
Requires  large  line  count  limit 

3 

Optional  output  3 . The  aerodynamic 
influence  coefficient  matrix,  the  right 
side  of  the  matrix  equation,  and  all 
solution  iterations 

4 

Optional  output  4.  This  option  prints 

out  the  successive  solution  iterations 
only 


N.B.:  The  normal  output  is  always  printed  in  addition  to  any  optional  out- 

put selected. 


Card  9A  - Single  Panel  Control  Card 


Column 

Variable 

Value 

Description 

1-10 

SINGPA 

0 

No  single  panel  input;  omit  card  set  9B, 

1 

continue  reading  input  cards.  Corner 
point  coordinates  of  this  panel  follow 
on  card  set  9B 

11-20 

NOPAN 

arbitrary 

Number  of  panels  to  be  deleted.  If  non- 

integer 

zero,  panel  indices  follow  on  card  set9C 
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Card  9B  - Panel  Corner  Point  Input 


Column 

Variable 

Value 

Description 

1-10 

X(I) 

arbitrary 

(floating 

point) 

x- coordinate 

of  corner  I 

11-20 

YU) 

II 

y-coordinate 

of  corner  I 

21-30 

Z(I) 

•1 

z- coordinate 

of  corner  I 

Repeat  card  9B  four  times,  once  for  each  corner  of  the  panel. 


Card  Set  9C  - Indices  of  Deleted  Panels  - NOPAN  indices  of  deleted  panels 
are  read  (7110  format)  if  NOPAN  > 0 on  card  9A.  A maximum  of  100  panels 
may  be  deleted.  Wing  and  body  vortex  lattice  control  panels  may  not  be 
deleted. 


Card  Set  10  - Body  Panel  Input  - This  card  set  allows  the  body  panels  to  be 
calculated  automatically  from  the  section  geometry  data.  Five  options  are 
available  for  inputting  the  section  geometry.  The  XYZ  program  input  refer- 
red to  below  conforms  with  the  format  of  Reference  2.  Omit  this  card  set 
if  CASE  = 2 on  card  8. 


Card  10A  - Number  of  Body  Sections 


Column 

Variable 

Value 

Description 

1-10 

NB 

arbitrary 

integer 

Number  of  body  sections  (2  NB  < 70) 

Card  10B 

- Body  Section  Geometry 

Column 

Variable 

Value 

Description 

1-10 

XBE 

arbitrary 

(floating 

point) 

x-coordinate  of  origin  of  body 
coordinate  system,  except  b3  ank 
program  input  format  is  used 

section 
when  XYZ 

11-20 

YBE 

II 

Simila-'ly  the  y-coordinate 

21-30 

ZBE 

II 

Similarly  the  z-coordj.nate 
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Column 

Variable 

Value 

Description 

31-40 

MB 

arbitrary 

integer 

Number  of  input  points  on  section 
(3  < MB  < 60, 

I • 

If  MB  < 0,  XYZ  program  input  format  re- 
quested 

50 

OPT 

0 

Body  section  geometry  input  by  y-z  co- 
ordinates on  card  set  10C  and  10D 

1 

Body  section  beometry  same  as  preceding 
body  section  - card  set  10C  or  10d  omit- 
ted. Note:  YBE  and  ZBE  are  additive  to 

preceding  values 

2 

Body  section  geometry  input  in  polar 
coordinates,  r,  9,  on  card  set  10C 

3 

Body  of  revolution,  section  geometry 
input  as  section  radius  and  theta  incre- 
ment on  card  set  10C 

60 

FLAG 

0 

Normal  body  section 

1 

Terminal  body  section  (end  of  current 
body  panel  network) 

2 

End  of  rectangular  grid 

Card  Set 

10C  - Body  Section  Coordinates  (Normal  Input) 

Column 

Variable 

Value 

Description 

1-10 

B ( J) 

arbitrary 
(f loating 
point) 

y-coordinate  of  point  J if  OPT  = 0;  or 
angular  coordinate  (in  degrees)  of  J 
if  OPT  = 2;  or  increment  angle,  A0,  in 

degrees  if  OPT  = 3 j 

11-20 

A ( J) 

z-coordinate  of  point  J if  OPT  =0;  or 
r-coordinate  of  point  J if  OPT  = 2;  or 
body  section  radius  if  OPT  = 3 

1 

21-30 

D(J) 

•• 

A x shift  of  point  J if  OPT  = 0 or  2 

Card  set  IOC  contains  MB  cards  if  OPT  = 0 or  2;  contains  only  1 
OPT  = 3;  and  is  omitted  if  OPT  = 1 or  MB  < 0 on  card  10B. 


card  if 
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Card  Set  10D  - Alternate  ZYZ  Input 


i 


Column 

Variable 

Value 

Description 

1-12 

D (J) 

arbitrary 

(floating 

point) 

x- coordinate 

of  point  J 

13-24 

B (J) 

If 

y- coordinate 

of  point  J 

25-  36 

A (J) 

z- coordinate 

of  point  J 

This  card  set  is  omitted  unless  MB  < 0 in  card  10B. 

Note:  Repeat  card  10B  and  card  sets  IOC  or  10D  NB  times  to  complete  card 

set  10. 


Card  llfl  - Number  of  Wing  Sections 
Column  Variable  Value 


1-10 


20 


MW 


KOORD 


arbitrary 

integer 


Description 

Number  of  wing  sections  (2  < NW  < 40) 


Wing  section  ordinates  input  in  percent 
of  local  chord 

Wing  section  ordinates  input  are  not 
normalized 


Card  11B  - Wing  Section  Geometry 


Column 

Variable 

Value 

Description 

1-10 

XBE 

arbitrary 

(floating 

point) 

x-coordinate  of  origin  of  wing  section 
coordinate  system 

11-20 

YBE 

II 

Similarly  the  y-coordinate 

21-30 

ZBE 

ft 

Similarly  the  z-coordinate 

31-40 

CHRD 

ft 

Chord  length  section 
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Cara’  11B  - Wing  Section  Geometry  (continued) 


Column 

Variable 

Value 

Description 

41-50 

ALF 

arbitrary 

(floating 

point) 

Section  twist  angle  (degrees)  or  flap 
rotation  angle.  (Degrees  positive  for 
positive  flap  deflection) 

51-60 

XAL 

II 

Center  of  twist  or  rotation  in  percent 
chord 

61-65 

MW 

arbitrary 

integer 

Number  of  coordinates  in  section 
(5  < MW  < 59) 

70 

OPT 

0 

Wing  section  ordinates  to  be  used  from 
card  set  11D 

1 

Wing  section  ordinates  same  as  preceding 
section  - card  set  11D  omitted 

75 

FLAG 

0 

Normal  case  - surface  vorticity  calcula- 
ted automatically 

1 

Terminal  wing  section  (end  of  current 
wing  panel  network) 

2 

No  vortex  lattice  panels  calculated  for 
this  section 

3 

The  coordinates  of  the  last  bound  vortex 
in  the  vortex  lattice  are  read  in  on 
card  13  for  this  section 

80 

DEL 

0 

No  wing  dihedral 

1 

Dihedral  input  on  card  set  11C 

Card  Set 

11C  - Wing 

Dihedral  Input 

Column 

Variable 

Value 

Description 

1-10 

DELTA 

arbitrary 

(floating 

point) 

Dihedral  angle  (degrees) 
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Card  Set  11C  - Wing  Dihedral  Input  (continued) 


Column 

Variable 

Value 

Description 

11-20 

YO 

arbitrary 

(floating 

point) 

y and  z ordinates  of  axis  of  rotation 
of  wing  panel 

21-30 

ZO 

II 

N.B.:  Omit  card  set  11C  if  DEL  = 0 in  card  11B. 


Card  Set  IIP  - Wing  Section  Coordinates 
Column  Variable  Value 


1-10 


11-20 


B ( J) 


A(J) 


arbitrary 

(floating 

point) 


Description 

x- coordinate  of  point  J 
2-coordinate  of  point  J 


N.B.:  Card  set  11D  contains  MW  cards  of  OPT  = 0,  and  is  omitted  if  OPT  = 1. 
Repeat  card  11B  and  card  sets  11C  and  11D  NW  times  to  complete  card  set  11. 

Card  Set  12  - Vortex  Lattice  Control  Point  Location 
Column  Variable  Value 


1-10 


11-20 


WAKE 


POINT 


arbitrary 

(floating 

point) 


Description 

Extension  of  vortex  lattice  into  wake  in 
percent  chord  (usually  100.) 


Location  of  vortex  lattice  control  point 
in  percent  chord  behind  trailing-edge 
(usually  . 1) 


N.B.:  These  values  are  not  used  if  FLAG  = 3 on  card  11B. 


Card  Set  13  - Relocation  of  Vortex  Lattice  Terminal  Points  - This  card  set 
is  omitted  unless  FLAG  = 3 on  card  11S.  For  each  wing  section  having  FLAG 
= 3,  two  additional  cards  are  required  to  specify  the  terminal  points  of 
the  streamwise  vortices. 
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Card  jet  13A  - Inboard  Terminal  Points 


Column 

Variable 

Value 

Description 

1-10 

XLP 

arbitrary 

(floating 

point) 

x- coordinate  of 
terminal  point 

inboard  edge 

of 

lattice 

11-20 

YLP 

If 

y- coordinate  of 
terminal  point 

inboard  edge 

of 

lattice 

21-30 

ZLP 

It 

z-coordinate  of 
terminal  point 

inboard  edge 

of 

lattice 

Card  Set  13B  - Outboard  Terminal  Points  - Same  as  card  13A  for  outboard  edge 
of  lattice  terminal  point. 


Card  Set  14  - Body  Vortex  Lattice  Input  - This  card  set  allows  additional 
vortex  lattices  to  be  located  inside  the  body  of  wing-body  combinations, 
and  is  omitted  if  CASE  < 3 on  card  8. 


Card  Set  14A  - Number  of  Streamwise  Vortices  in  Body  Vortex-Lattice  Network 
Column  Vai  able  Value  Description 


1-10 


NV  arbitrary  Number  of  streamwise  vortices  in  body 

integer  vortex- lattice  network  (NV  < 40) 


Note:  The  sum  of  all  wing  and  body  vortex- lattices  may  not  exceed  35. 


Card  Set  14B  - Vcrcpx- Lattice  Geometry 


Column 

Variable 

Value 

Description 

1-10 

XBE 

arbitrary 

(floating 

point) 

x-coordinate 

vortex 

of  origin  of 

streamwise 

11-20 

YBE 

II 

y- coordinate 
vortex 

of  origin  of 

streamwise 

21-30 

ZBE 

It 

z-coordinate 

of  origin  of 

streamwise 

vortex 
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31-40 

MV 

arbitrary 

integer 

Number  of  bound  vortices  in  lattice 
2 < MV  < 60 

41-50 

OPT 

0 

Vortex- lattice  points  to  be  read  from 
card  set  14C 

1 

Vortex  lattice  points  same  as  preceding. 
Omit  card  set  14C 

2 

Optional  vortex-lattice  control  panel 
coordinates  read  on  card  14D3 

51-60 

FLAG 

0 

Normal  case  - vortex-lattice  panels  cal- 
culated 

1 

Terminal  vortex  of  current  body  vortex- 
lattice  network 

2 

Corner  points  of  control  point  panel  to 
be  read  on  cards  14D2  and  14D3  (used 
when  arbitrary  control  point  is  desired) 

61-70 

SIMPOT 

0 

Symmetry  option  specified  on  card  8 
enforced  for  this  vortex 

1 

Symmetry  option  ignored  for  this  vortex 
lattice  (used  for  inserting  vortex- 
lattice  networks  in  vertical  tails  loca- 
ted in  x,z  plane) 

Card  Set 

14C  - Vortex- 

-Lattice  Coordinates 

Column 

Variable 

Value 

Description 

1-10 

B ( J) 

arbitrary 

(floating 

point) 

x-coordinate  of  point  J 

11-20 

A ( J) 

II 

z-coordinate  of  point  J 

21-30 

C(J) 

•1 

Vortex-lattice  strength  at  point  J 
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Card  Set  X4C  - Vortex- Lattice  Coordinates  (continued) 


Column  Variable  Value 


Description 


31-40  D(J)  arbitrary  Ay  shift  of  point  J 

(floating 

point) 


Card  set  14C  contains  MV  cards  if  OPT  = 0,  and  is  omitted  if  OPT  = 1 on 
card  14B. 


Control  Set  14D  - Vortex-Lattice  Terminal  Point  and  Control  Point  Coordinates 
Two  or  three  additional  cards  are  required  to  specify  the  terminal  point  of 
the  streamwise  vortex,  and  the  corner  points  of  the  lattice  control  point 
panel . 


Card  14D1 


Column 

1-10 


11-20 


Variable  Value 


Description 


B arbitrary  x-coordinate  of  terminal  point  of  stream- 

(f loafing  wise  vortex 

point) 


A 


z-coordinate  of  terminal  point  of  stream- 
wise  vortex 


21-30  D 


Ay  shift  of  terminal  point  of  streamwise 
vortex 


Note:  This  point  also  defines  the  upstream  corner  of  the  control  point 

panel  if  FLAG  / 2 on  card  14B. 


Card  14D2  - Same  as  card  14D1,  containing  the  coordinates  of  the  downstream 
comer  of  the  control  point  panel  if  FLAG  ? 2 on  card  14B.  If  FLAG  = 2 
on  card  14B,  this  card  contains  the  coordinates  of  the  upstream  corner  of 
the  control  point  panel. 

Card  14D3  - If  FLAG  = 2 on  card  14B,  this  card  contains  the  coordinates  of 
the  downstream  comer  of  the  control  point  panel  in  the  same  format  as  card 
14D1.  Omit  this  card  if  FLAG  / 2 in  card  14B1. 

Note:  Repeat  card  14B  and  card  sets  14C  and  14D  NV  times  to  complete  card 

set  14. 


84 


Aerodynamic  Input  Cards  - The  configuration  panel  geometry  is  transferred 
to  the  aerodynamic  section  of  the  program  by  TAPE  11.  Additional  aerody- 
namic input  cards  required  are  described  below. 


Card  15  - Case  Identification  Card  - Card  15  contains  any  desired  case  iden- 
tification in  columns  1-80. 


Card  16  - Iteration  Option  Card 

Column  Variable  Value 

1-10  NIT  arbitrary 

integer 

11-20  I EPS 

21-30  I TYPE  1 

2 


Description 

Maximum  number  of  iterations  (15  - 25) 


Exponent  of  10  setting  limiting  value 
for  residue  of  iterative  solution  (-3  or 
-4  recommended) 

Blocked  Jacobi  iteration  procedure 
Blocked  Gauss-Siedel  iteration  procedure 


Blocked  Gauss-Seidel  with  controlled 
successive  over- re taxation 


4 Blocked  Gauss-Seidel  with  successive 

over- relaxation 


N.B.:  Options  3 and  4 have  not  been  completely  checked  out,  but  will  be 

made  available  at  a later  date. 


Card  17 

- Configurat 

ion  Options 

Column 

Variable 

Value 

./ 

Description 

1-10 

COMPT 

0 

Forces  and  moments  calculated  for  com- 
plete configuration 

arbitrary 

integer 

Forces  and  moments  calculated  on  com- 
ponents. Panel  indices  of  each  compo- 
nent follow  on  card  26 

11-20 

SECT 

0 

No  wing  section  forces  and  moments 
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Card  17  - Configuration  Options  (continued) 


Column 

Variable 

Value 

Description 

11-20 

SECT 

1 

Wing  section  forces  and  moments  calcula- 
ted. Wing  section  indices  follow  on 
card  25,  panel  indices  in  each  section 
on  card  26,  and  section  reference  lengths 
on  card  29 

2 

Forces  and  moments  calculated  on  subsec- 
tions. The  number  of  subsections  follow 
on  card  25,  the  number  of  panel  groups 
on  card  27,  the  panel  indices  in  each 
group  on  card  28,  and  subsection  refer- 
ence lengths  on  card  29 

Card  18  - 

Reference 

Parameters 

Column 

Variable 

Value 

Description 

1-10 

REFA 

arbitrary 

(fixed 

point) 

Reference  data 

11-20 

REFL 

It 

Reference  chord  (MAC) 

21-30 

XOO 

II 

Axial  distance  of  leading-edge  of  MAC 
from  origin 

31-40 

X25 

II 

Axial  distance  of  quarter  chord  of  MAC 
from  origin 

N.B.:  Default  option  — all  variables  set  to  1.  internally  (leave  card 

blank) . 


Card  19  - 

Configuration 

Lift 

Option 

Column 

Variable 

Value 

Description 

1-10 

KUT 

0 

Nonlifting  configuration,  no  vortex-lat 
tice  Kutta  condition  imposed 

1 

Lifting  configuration,  vortex-lattice 
Kutta  condition  imposed 
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Card  19  - Configuration  Lift  Option  (continued) 


Column 

Variable 

Value 

Description 

1-10 

KUT 

-1 

Wing  vortex  lattice  extends  through  body 
having  same  strength  as  adjacent  wing 
vortex  lattice 

11-20 

NBV 

arbitrary 

integer 

Number  of  body  vortices  (NBV  < 5) 

21-30 

NV  (1) 

" 

Number  of  wing  vortices  associated  with 
body  vortex  1 

31-40 

NV  (2) 

" 

Number  of  wing  vortices  associated  with 
body  vortex  2 

61-70  NV(5)  " Number  of  wing  vortices  associated  with 

body  vortex  5 

Card  20  - Compressibility  Rule  Option 


Column 

Variable 

Value 

Description 

1-10 

KOMPR 

1 

Gothert 

Rule  1 selected 

2 

Gothert 

Rule  2 selected 

11-20 

POINTS 

0 

One  body 
points 

component  case,  no  off-oody 

1 

Flap  case,  off-body  points  will  be  cal- 
culated 

Card  21  - 

Number  of 

Mach  Numbers 

Column 

Variable 

Value 

1 

Description 

1-10 

NMA 

arbitrary 

integer 

Number  of  Mach  numbers  following  on  card 

set  22 
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Card  22 

- Mach  Number 

Column 

Variable 

Value 

Description 

1-10 

MA 

arbitrary 

(floating 

point) 

Mach  number 

Card  2 3 

- Number  of 

Angles-of-Attack  and  Yaw 

Column 

Variable 

Value 

Description 

1-10 

NAL 

arbitrary 

integer 

Number  of  angles-of- attack  following  on 
card  set  24 

Card  24 

- Anqle-of  Attack  and  Yaw 

Column 

Variable 

Value 

Description 

1-10 

ALPHA 

arbitrary 

(floating 

point) 

Angle-of- attack  in  degrees 

11-20 

BETA 

II 

Angle  of  yaw  in  degrees 

Card  25 

- Number  of 

Sections 

Column 

Variable 

Value 

Description 

1-10 

IS 

arbitrary 

integer 

Number  of  sections;  omit  if  SECT  = 0 on 
card  16 

Card  26 

- Panel  Indices 

Column 

Variable 

Value 

Description 

1-10 

IA 

arbitrary 

integer 

Index  of  initial  panel  in  section 

11-20 

IE 

II 

Index  of  final  panel  in  section 
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Card  27  - Number  of  Panels  in  Subsections 


Column 

Variable 

Value 

Description 

1-10 

IREI 

arbitrary 

integer 

Number  of  panels  in  subsections.  Omit 
if  SECT  < 2 on  card  17 

Card  28  - 

Subsection 

Panel  Indices 

Column 

Variable 

Value 

Description 

1-5 

11(1) 

arbitrary 

integer 

Panel  indices  of  all  panels  in  subsec- 
tion; omit  if  SECT  < 2 on  card  17 

6-10 

11(2) 

11 

11-15 

11(3)  . 

..  etc. 

Card  29  - 

Reference 

Lengths 

Column 

Variable 

Value 

Description 

1-10 

DELY 

arbitrary 

integer 

Width  of  section 

11-20 

REFL 

II 

Reference  length  of  section 

21-30 

XLE 

•1 

Moment  reference  point  of  section 

N.B.:  Cards  25  - 29  must  be  repeated  for  each  angle-of-attack  or  yaw,  if 

section  data  requested. 

Streamline  Input  Cards  - Each  streamline  selected  is  constrained  to  pass 
through  a chosen  panel.  This  avoids  the  possibility  of  the  selected 
streamlines  bunching  together  in  the  separated  flow  region. 

Card  30  - Case  Identification  Card  - Card  30  contains  identification  for 
the  streamline  calculations. 


89 


Card  31  - Streamline  Option  Card 


Column  Variable  Value  Description 

1-10  ENDPILE  -1  Program  switches  out  of  WBAERO  onto  the 

streamline  overlay 

Card  32  - Streamline  Panel  Selection 


Column 

Variable 

Value 

Description 

1-5 

NLIN 

variable 

Number  of  prescribed  panels 
line  calculations 

for  stream- 

6-10 

NSP(l) 

variable 

Panel  number  for  streamline 

number 

1 

11-15 

NSP (2) 

»» 

Panel  number  for  streamline 

number 

2 

16-20 

NSP (3) 

. . . etc. 

N.B.  : 

Repeat  this 

card  TRMAX 

times . 

Card  33  - End  of  Data  - 6/7/8/9 
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PROGRAM  OUTPUT  DESCRIPTION 


The  standard  output  of  the  program  is  a function  of  iteration  number  as 
described  in  the  following  paragraphs. 

In  Iteration  1,  tables  of  parameter  values  for  variably  dimensioned  arrays 
from  subroutines  WBOLAY  and  WBPAN  are  printed  out,  followed  by  a table  of 
panel  corner  points  for  the  input  geometry.  Output  from  subroutine  WBAERO 
is  printed  next;  the  list  of  variable  dimensioned  array  parameters  for  this 
subroutine  is  followed  by  selected  input  parameters  and  a list  of  central 
processor  times.  The  first  time  is  at  the  entry  point  to  the  subroutine* 
and  the  difference  between  the  second  and  third  times  gives  the  time  to 
form  the  matrix  of  influence  coefficients.  The  fourth  time  is  printed  just 
before  entering  the  SOLVE  sub routine, and  the  fifth  time  (after  SOLVE)  is 
printed  out  after  SOLVE  has  given  the  history  of  the  residual  in  the  itera- 
tive solution  method  and  the  solved  singularity  strengths  (in  order  of  panel 
number) . 

The  next  table  gives  the  attached  flow  velocity  vector  and  magnitude  and 
the  pressure  coefficient  at  each  panel  control  point,  and  is  followed  by  a 
summary  of  the  total  coefficients  which  include:  the  axial,  normal  and 

side  forces;  the  moments  about  the  coordinate  axes;  the  pitching  moment 
about  the  reference  point  and  about  the  quarter  chord  of  the  MAC;  the  x-wise 
center  of  pressure;  and  the  lift  side  force  and  drag  coefficients. 

All  the  following  output  is  repeated  for  each  iteration. 

Output  from  the  streamline  program  starts  with  the  variable  dimensioned 
array  parameters.  The  prescribed  streamline  starting  points  and  panels  are 
given  next  followed  by  the  tabulated  characteristics  for  each  streamline. 

The  tables  list  the  x,y,z  coordinates,  the  velocity  vector  and  magnitude, 

the  pressure  coefficient,  the  geodesic  curvature  (K  , K ) , the  metric  coef- 

1 2 

ficient  (H  ) and  the  distance  from  the  upstream  end  of  the  calculated 
2 

streamline . 

Next,  the  calculated  boundary  layer  characteristics  are  tabulated  for  each 
streamline  and  include  the  local  shape  parameter  and  skin  friction  coef- 
ficient. Each  table  is  preceded  by  a summary  of  the  streamline  geometry 
and  pressure  distribution. 

Subroutine  WAKECON  output  follows  the  boundary  layer  characteristics  and 
starts  with  a summary  of  the  separation  points  and  the  separation  panels 
that  give  an  approximate  representation  of  the  separation  line.  The  geo- 
metric parameters  for  the  separated-wake  vorticity  panels  are  then  tabulated, 
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and  are  followed  by  the  singularity  strength  solution:  surface  sources 

(sigma)  and  wake  vorticity  (gamma) . 

Tables  of  velocity  and  pressure  coefficients  in  the  presence  of  the  separa- 
tion model  are  the  last  output  for  each  iteration.  This  output  starts  with 
the  mean  velocity  (UM) , the  pressure  coefficient  (CP) , and  the  wake  total 
head  decrement  (DH)  at  the  wake  panel  control  points,  and  is  followed  by 
the  surface  velocity  vector  and  magnitude  and  the  pressure  coefficient  at 
each  panel  control  point. 

After  the  last  iteration,  the  wetted  area,  the  reference  area,  and  the  drag 
summary  table  are  printed  out;  the  table  includes  the  drag  coefficients  for 
the  skin  friction,  the  pressure,  and  the  total  based  on  the  reference  area. 
The  total  drag  is  also  presented  in  D/q  form  in  the  units  of  the  reference 
area. 
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LIST  OF  SYMBOLS 


B.  . 

ID 


Aerodynamic  influence  coefficient 
Normal  velocity  due  to  external  source 
Total  drag  coefficient 

Profile  drag  coefficient  (Squire  and  Young) 
Base  pressure  drag  due  to  separation 


C Lift  coefficient 

Xj 

C„  Moment  coefficient 

M 

Local  skin  friction  coefficient 

C Skin  friction  drag 

Df 

C Pressure  coefficient 

P 

C Shear  stress  integral 

c Chord 

F,  F , G,  G Universal  functions  in  Curie's  laminar  method 
o o 

g Correction  term  to  Thwaites'  laminar  method 

H Shape  factor,  ratio  of  displacement  to  momentum  thickness 

(6*/9) 

K Non-dimensional  pressure  gradient  parameter 

M Mach  number 

M Local  Mach  number 

L 
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LIST  OF  SYMBOLS  (Continued) 


P 

R 


®e 

^ins 

R0trans 


r 

S 

u 


u 


T 


V 


X#  y r Z 


a 

e 

Y 


Free  stream  Mach  number 

Number  of  singularities 

Total  normal  velocity  at  it*1  point 

Static  pressure,  pounds  per  square  inch  absolute 
Local  radius  of  curvature 
Chord  Reynolds  number  U^c/V 

Momentum  thickness  Reynolds  number  U0/V 

Streamwise  momentum  thickness  Reynolds  number  at  insta- 
bility point 

Streamwise  momentum  thickness  Reynolds  number  at 
transition 

Transverse  radius  of  curvature,  inches 
Distance  along  a streamline 
Local  velocity  at  edge  of  boundary  layer 
Free  stream  velocity 

V2 

Friction  velocity  (x  /p) 

w 

Tangential  velocity  at  body  surface 
Cartesian  coordinates  of  point 
Angle-of- attack 
Angle  of  yaw 
Vortex  strength 
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LIST  OF  SYMBOLS  (Continued) 


6 Boundary  layer  thickness 

p Density  of  air 

0 Source  strength 

T Shear  stress 

T Local  surface  shear  stress 

w 


Subscripts 

i 

. th 

i value 

in 

Incompressible 

ins 

Instability 

j 

.th 

j value 

L 

Local  Value 

1 

Lower 

trans 

Transition 

t 

Turbulent 

u 

Upper 

i 
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